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Abstract 


The focus of this work is the buckling response of symmetrically laminated composite plates having a 
planform area in the shape of an isosceles trapezoid. The loading is assumed to be inplane and applied 
perpendicular to the parallel ends of the plate. The tapered edges of the plate are assumed to have simply 
supported boundary conditions, while the parallel ends are assumed to have either simply supported or 
clamped boundary conditions. A semi-analytic closed-form solution based on energy principles and the 
TrefFtz stability criterion is derived and solutions are obtained using the Rayleigh-Ritz method. Intrinsic in 
this solution is a simplified prebuckling analysis which approximates the inplane force resultant distributions 
by the forms Nx=P/W(x) and Ny=Nxy=0, where P is the applied load and W(x) is the plate width which, for 
the trapezoidal planform, varies linearly with the lengthwise coordinate x. The out-of plane displacement is 
approximated by a double trigonometric series. This analysis is posed in terms of four nondimensional 
parameters representing orthotropic and anisotropic material properties, and two nondimensional parameters 
representing geometric properties. For comparison purposes, a number of specific plate geometry, ply 
orientation, and stacking sequence combinations are investigated using the general purpose finite element 
code ABAQUS. Comparison of buckling coefficients calculated using the semi-analytical model and the 
finite element model show agreement within 5%, in general, and within 15% for the worst cases. In order to 
verify both the finite element and semi-analytical analyses, buckling loads are measured for graphite/epoxy 
plates having a wide range of plate geometries and stacking sequences. Test fixtures, instrumentation 
system, and experimental technique are described. Experimental results for the buckling load, the buckled 
mode shape, and the prebuckling plate stiffness are presented and show good agreement with the analytical 
results regarding the buckling load and the prebuckling plate stiffness. However, the experimental results 
show that for some cases the analysis underpredicts the number of halfwaves in the buckled mode shape. In 
the context of the definitions of taper ratio and aspect ratio used in this study, it is concluded that the 
buckling load always increases as taper ratio increases for a given aspect ratio for plates having simply 
supported boundary conditions on the parallel ends. There are combinations of plate geometry and ply 
stacking sequences, however, that reverse this Irend for plates having clamped boundary conditions on the 
parallel ends such that an increase in the taper ratio causes a decrease in the buckling load. The clamped 
boundary conditions on the parallel ends of the plate are shown to increase the buckling load compared to 
simply supported boundary conditions. Also, anisotropy (the D16 and D26 terms) is shown to decrease the 
buckling load and skew the buckled mode shape for both the simply supported and clamped boundary 
conditions. 
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Chapter 1 
Introduction 


1.1 Description of Problem 

Thin plates are fundamental components of many engineering structures. The aerospace industry exploits 
the high stiffness-to-weight and strength-to-weight ratios afforded by fiber reinforced composite plates by 
using them as both primary and secondary structures to reduce the weight of high performance aircraft. 
Many of these plates, particularly wing skins and fuselage panels, are used in stiffness critical applications. 
Thus, for those plates loaded in compression, buckling must be considered as the primary mode of failure. 
Many of the composite plates on aircraft, again notably wing skins, have a nonrectangular planform area; in 
fact, many of the composite plates in aircraft applications have a planform area that is tapered along the 
length. The purpose of this study is to investigate the buckling response of composite plates having planform 
area in the shape of an isosceles trapezoid. 

Figure 1 . 1 presents the planform geometry of an isosceles trapezoid considered in this study. The plate is 
defined to have length L, widths W! and W 2 , and half-widths b, and t >2 on the narrow and wide ends, 
respectively. The trapezoid is considered isosceles because the comer angle 0 is the same for both the right 
and left halves of the plate. The coordinate system originates in the center of the narrow end of the plate. The 
x axis coincides with the centerline of the plate, and the y axis traverses the narrow end. 
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Figure 1.1 Trapezoidal Plate Geometry 


1.2 Survey of Literature 

An investigation of the buckling response of plates having a tapered planform first appeared in the open 
literature in 1956. Klein [1] presented design curves for determining the buckling loads of simply supported 
isotropic plates having the shape of an isosceles trapezoid. Simply supported boundary conditions were 
assumed on all four edges, axial compression loads were imposed on the parallel ends, and shear loads were 
assumed to act along the tapered edges to allow the plate to remain in static equilibrium for any ratio of axial 
loads on the ends of the plate. The analysis used the method of collocation in which a deflected shape is 
assumed, appropriate derivatives are calculated, the derivatives are evaluated numerically at several discrete 
collocation points, and these equations are then substituted into the governing differential equations to yield 
an eigenvalue problem. The method of collocation does not required the assumed displacement functions to 
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be integrable in closed form, form an orthogonal set, or even satisfy the boundary conditions exactly for the 
entire plate. The method only requires the assumed shape to be reasonable at the collocation points. Klein 
used a trigonometric function for the assumed deflected shape and developed a solution using three 
collocation points distributed along the centerline of the plate. The results showed an increase in the 
buckling load of a plate with a tapered planform as compared to the buckling load of a rectangular plate of 
the same length and with the same width as the wide end of the trapezoidal plate. No experimental results 
were included for verification. 

In 1956 Klein [2] also investigated the buckling response of plates tapered in both planform and thickness. 
The method of collocation was again used for the analysis. Many graphs were presented to show the effect 
of plate geometry on buckling over a very wide range. Again, however, experimental verification was not 
included. 

In 1957 Klein [3] presented a method of solution for the shear buckling of simply supported isotropic plates 
tapered in planform. The analysis was again based on the method of collocation using three points along the 
centerline of the plate. The buckling load and deflected mode shape were presented for one specific 
geometry. No experimental results were included for verification. 

Pope [4], in 1962, studied the effects of different boundary conditions on the buckling loads of isotropic 
plates having a planform area in the shape of an isosceles trapezoid. This analysis was developed using the 
Rayleigh-Ritz method and assumed that the buckled shape across the width of a tapered plate differs little 
from the buckled shape across a rectangular plate under uniform end load and with the same boundary 
conditions along the sides. Results are presented graphically for plates with opposite pairs of edges either 
simply supported or clamped. These results showed an increase in the buckling load when the planform is 
tapered as compared to the buckling load of a rectangular plate of the same length and with the same width 
as the wide end of the trapezoidal plate. Clamping the ends of the plate is also shown to increase the 
buckling load. For comparison purposes, the buckling loads presented by Pope are usually lower than those 
presented by Klein. It is somewhat disturbing, however, that Pope doesn’t obtain the classical solution for 
the buckling loads for rectangular uniaxially loaded isotropic plates having simple supports on all four edges 
[5], a festooned curve with the minimum value of 4 at aspect ratios 1, 2, 3 etc. (if the buckling parameter has 
the classical definition of N x b 2 /7t 2 D). The results presented by Pope have a minimum value of 
approximately 3.7 at aspect ratios 1, 2, 3, and 4. 

The buckling analyses for plates having tapered planform presented in the technical literature to date are 
restricted to isotropic plates. There is a significant amount of information available regarding the buckling 
response of composite plates, but this is limited to rectangular geometries. Of particular importance in a 
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buckling analysis when considering composite plates is the effect of the anisotropic plate bending stiffness 
terms D 16 and D 2 5 . These terms couple the out-of-plane bending and twisting of the plate. Because these 
stiffness terms appear as coefficients of terms involving mixed partial derivatives of the displacement 
functions, including these terms greatly increases the complexity of an analysis and hence they are often 
neglected. Neglecting these terms, however, can yield significant errors for certain anisotropic plates 

In 1969 Chamis [6] presented a criterion, based on experimental results, for determining when the 
anisotropic bending stiffness can be neglected. Nemeth [7], in 1986, presented nondimensional parameters 
that can be used to assess when the anisotropic bending stiffness can be neglected. These nondimensional 
parameters resulted naturally from the nondimensionalization of the buckling equation. Finite element 
results showed that a specially orthotropic analysis (D 16 =D 26 = 0 ) can yield nonconservative errors as high as 
25% in the predicted buckling loads for certain anisotropic cases. 

Nemeth [8] used the same nondimensional formulation of the buckling equation in 1992 to study long 
symmetrically laminated plates subjected to compression, shear, and inplane bending loads. Variational 
methods were used to derive the buckling analysis in terms of nondimensional parameters that characterize 
the plate bending orthotropy and plate bending anisotropy. Because the analysis is constructed in this 
manner, it is well suited for parametric studies. Results are presented as families of parametric curves that 
cover a wide range of material properties, plate geometries, and load combinations. 

1.3 Objective and Scope 

In general, the survey of literature shows that the buckling of plates tapered in planform is limited to 
isotropic cases, while the buckling of composite plates is limited to rectangular geometries, often with the 
influence of D 16 and D 26 neglected. Therefore, the objective of this study is to determine the buckling 
response of composite plates which have a trapezoidal planform area and which are subjected to uniform 
end-shortening. The influence of D 16 and D 26 are included. The investigation is conducted using three tools. 

• A special purpose analysis is developed using variational energy methods and the Rayleigh-Ritz 
method of solution. The analysis is posed in nondimensional parametric form, well suited for 
isolating and investigating the effects of individual material and geometric properties. 

• A finite element analysis is performed to verify several assumptions made in developing the 
Rayleigh-Ritz analysis, and to compare results for some specific plate geometries. 

• Buckling loads are determined experimentally for a series of graphite/epoxy composite specimens 
to verify both the finite element and Rayleigh-Ritz analyses. 
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The scope of this study is limited to symmetrically laminated composite plates having simply supported 
boundary conditions on the nonparallel edges, while the parallel edges have either simply supported or 
clamped boundary conditions. Known inplane compressive loads applied perpendicular to the parallel edges 
are assumed to represent the effects of uniform end shortening. 

The question may be asked: “Why bother with a semi-analytic approach if solutions can be obtained using 
the finite element method?” Because there is a wide range of materials, ply orientation angles, stacking 
sequences, and plate geometries, it is useful to present buckling loads in the form of design charts consisting 
of nondimensional parametric relationships. These relationships can present buckling coefficients over a 
wide range of materials, geometries, and boundary conditions in relatively few graphs. Additionally, these 
parametric relationships can be grouped into sets to reveal trends in the structural response governed by 
certain parameters. Because results must be calculated over a wide range of combinations of design 
parameters when generating these design charts, a special purpose analysis is preferred over a general 
purpose finite element analysis due to the cost and effort involved in developing numerous finite element 
models. Also, on a more philosophical level, developing the analysis using energy methods allows one to 
“keep their fingers in the physics” and investigate the effects of various parameters at many stages 
throughout the analysis before buckling loads are ever calculated. 

In what follows, Chapter 2 presents the basic assumptions, develops the analysis based on variational energy 
methods, presents the assumed shape for the out-of-plane displacement, and discusses the computer 
implementation. Chapter 3 discusses the case of simply supported boundary conditions on the parallel ends 
of the plate and presents figures comparing analytical results with finite element predictions. A series of 
parametric curves is then presented for this set of boundary conditions. Chapter 4 discusses the case of 
clamped boundary conditions on the parallel ends of the plate and presents figures comparing analytical 
results with finite element predictions. A second series of parametric curves is then presented for this set of 
boundary conditions. Chapter 5 describes the experimental method used to measure the buckling load and 
Chapter 6 presents the experimental results for comparison with the analytical results. The appendices 
contain much of the details of the mathematics involved in Chapter 2. 
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Chapter 2 

Development of Analysis 


2.1 Overview of Method 

A nondimensional analysis for the buckling loads of symmetrically laminated plates having trapezoidal 
planform and subjected to inplane uniaxial compression is developed in this chapter. The analysis is based 
on the Rayleigh-Ritz method. This semi-analytic method assumes a series approximation for the out-of- 
plane, or buckling, displacement of the plate. The analysis originates by formulating the total potential 
energy of a plate subjected to inplane loading. Because the Trefftz criterion argues that a change in the 
stability of a structure occurs when the first variation of the second variation of total potential energy of the 
structure equals zero, the first variation of the second variation of total potential energy is computed and then 
expressed in a nondimensional form suited for parametric studies. A series approximation satisfying the 
kinematic, i.e., displacement, boundary conditions of the plate is substituted into the nondimensional 
expression of the first variation of the second variation of the total potential energy. Assuming that the 
resulting surface integrals can be performed on a term-by-term basis yields a symmetric general eigenvalue 
problem that can be solved using standard routines on a computer. The prebuckling equilibrium conditions 
are addressed directly by assuming an inplane force resultant distribution. The prebuckling equilibrium 
equations are not solved explicitly. 

2.2 Basic Assumptions 

The development of the buckling analysis is based on the following basic assumptions: 

(1 ) Each layer in the laminate is orthotropic, linear elastic, and of constant thickness. 

(2) The plate thickness is small compared to its length and width. 

(3) The plane-stress assumption is valid: the out-of-plane normal and shear stresses are zero. 

(4) Kirchhoff’s assumption is valid: line elements normal to the reference surface of the plate 
remain normal and inextensible during deformation. 
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(5) Von Karman nonlinear strain -displacement equations are applicable: strains are small 
compared to unity, out-of-plane rotations are moderate, and rotations about the normal to the 
reference surface are zero. 

(6) No body forces exist. 

(7) The plates are initially flat. 

Nonlinear strain-displacement equations are necessary because when an initially flat plate loses stability, the 
out-of-plane displacements are no longer linearly proportional to the applied inplane loading. The von 
Karman plate theory accounts for this by allowing moderate transverse rotations during displacement. This 
yields a geometrically nonlinear relationship between the inplane forces and out-of-plane displacements 
because the plate membrane forces are contributing to the out-of-plane equilibrium of the plate. 


2.3 Total Potential Energy 

The total potential energy of a plate is defined as 


H/2 


n(u°, V°, W°) = j { a x £ x + CT y £ y + a z e z + T xyY*y + T *^« + T y^ }dZdA 


( 2 . 1 ) 


A -H/2 


. rr / ° 0 °\ 

+ n ioad( U - V - W ) 


where n }oad is the potential energy due to the tractions on the surfaces of the plate, and the superscript o 
indicates quantities defined at the reference surface of the plate. With the plane-stress assumption, Eq. (2.1) 
can be simplified to 


H/2 


I — i / o o O \ 

n (u , v , w ) = - 


= \ J J J i CT x £ x + °y £ y + X xy Yxy > dzd A + n .oad 


U U U 

(u , V , w ) 


(2.2) 


A -H/2 


Using Kirchhoff’s assumption, the strains throughout the thickness of the plate can be expressed in terms of 
the reference surface strains and curvatures as 


£ 

£ 


X 

y 


Yxy 


= £ + ZK 


= e y + Z Ky 


= y + zk 
**y *y 


(2.3) 
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Substituting Eq. (2.3) into Eq. (2.1) yields 


H/2 


FT , U O O 

n (u , v , w ) 


-ujj* 


(£%Z<) +0 y( £ y + ZK y) + ( Y^y + } dzd A 


(2.4) 


A -H/2 


1 -r y 0 ° 0 . 

+ n ioad( U > V - W ) ' 

Integrating through the thickness and using the following standard definitions for the inplane force and out- 
of-plane bending moment resultants, 


H/2 


-H/2 

H/2 


N >- J 


Oydz 


-H/2 

H/2 


-H/2 


H/2 

N ,= | a * dz M * = | ° x zdz 


-H/2 

H/2 


M y = J o y zdz 


-H/2 

H/2 


N *y= J T *y dZ M *y=J\y Zdz ’ 


-H/2 


(2.5) 


yields 


° O 0 

n (u , v , w ) = - 


=ifj 


(N x e° + N y eJ + N iy Y° y + + M y K° + M xyK “ y ) dA 


(2.6) 


T-T . W U U 

+ n .oad(^ V ’ W ) 


2.4 Variation of Total Potential Energy 

As indicated by the notation of Eq. (2.6), the total potential energy of a plate is a function of the kinematic 
variables u°, v°, and w°. To develop the stability equations, the variation of the total potential energy with 
respect to arbitrary but infinitesimal variations in these kinematic variables must be investigated. The total 
potential energy of a plate computed with arbitrary infinitesimal variations in the displacements u°, v°, and 
w° can be written as 
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n + An = ^JJ{(N x + AN x )( e %Ae°) + (N y + AN y )(e; + Ae;) (2.7) 

A 

+ (N xy + AN xy ) (y° xy + Ay xy ) + (M x + AM x ) (k° + Ak°) 

+ (M y + AM y ) (K° y + AK° y ) + (M xy + AM xy ) (< y + AK° y ) }dA 
+ n + An ioad (u° + Au°, v° + Av°, w° + Aw°) ’ 

where the notation A ( ) represents an increment in ( ) due to the increments Au°, Av°, and Aw°. At this 
point it is necessary to express the right hand side of Eq. (2.7) completely in terms of the displacements u , 
v°, and w° and their increments. It is important to note that both the displacements and their increments 
satisfy the kinematic boundary conditions of the plate. The strains in Eq. (2.7) are calculated by substituting 
the displacements into the strain-displacement equations, and the stress resultants are computed by 
substituting the strains into the constitutive equations and integrating with respect to z, as specified by Eq. 
(2.5). 


Using the von Karman assumption of moderate rotations, the reference surface strains and curvatures are 


defined as 


and 



( 2 . 8 ) 


(2.9) 


Using these strain-displacement relationships, the (e x + AeJ term in Eq. (2.7), for example, can be 
expanded as follows: 
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£° + Ae° = — (u° + Au°) + - (w° + Aw°)~| . (2.10) 

x x Bx 2 LBx J 

It is more convenient to express the variations, or increments, in the displacement field using small 
parameter notation, namely, 


. ° O 

Au = eu, 

Av° = £v° (2.11) 

.0 0 

Aw = £W| • 

Using this notation, the displacements Up v ° are functions of x and y which must satisfy the same 

kinematic boundary conditions as u°, v°, and w°. The subscript 1 indicates that these are displacements in 
the neighborhood of the equilibrium displacement field, and the scalar parameter £ is assumed to carry the 
infinitesimal smallness of the variation. Consequently, the first, second, third, and fourth variations of the 
responses may be determined by grouping terms with like power of £, i.e., £, £ 2 , £ 3 , £ 4 . Using the notation of 
Eq. (2.11), Eq. (2.10) can be expressed as follows: 


Therefore, 


o .o r),o o. 1 Pi) / o 

£ * + 4£ > = ^ ( “ +E "' ,+ 2fc (W+£W ‘ ) ] 

2 ( \ 2 

Bx 2 ^ Bx ' ^ Bx 3x 9x ' 2 V Bx > 

o o 2 o 
= £ x + ££ x + £ £ - 


A£ 


o 

X 


o 2 o 

££ + £ £ , 
X, 


where 


o _ Bu? Bw^Bw^ 

X] Bx Bx Bx 

dx ) 



In a similar fashion, the (e° + A£°) 


term can be expanded as 


( 2 . 12 ) 


(2-13) 


(2.14) 
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e° + Ae° = ^-(v° + Av°) + irf(w° + Aw°)l 
y y dy 2 Ldy J 

2 

°v l r r) / ° 

= ^ (V+£V|)+ ife <w+ewi) ] 

dy 2 v dy ' ^ dy dy dy ' 2 

o o 2 o 

= e v + ee + e e . 
y y i y 2 


2 



( 2 . 15 ) 


Then, 


where 


A £„ 


o 2 o 

ee v + ££ v> 

*1 y 2 


o = dVj C^dw? 

y * 5y 3y dy 
2 

)• 

Finally, the shear strain term can be expanded as 


e°=ife 

y > 2 ^ 3 v 


7 ° + Ay° = — (u° + Au°) + — (v° + Av°) + — (w +Aw )^-(w + Aw ) 
xy xy dy dx dx dy 

= -^-(u° + EU°) + — (v° + EV°) + A( w ° + ewj)^(w° + ew°) 
dy dx dx dy 

du° dv° dw°dw° (du° 3v° dw°9w° 3w®3w°V 

9y dx dx dy ^ dy dx dx dy dx dy 

^ dx dy 

As a result. 


A V*y = £Y xy] 


2 

+ e v 


*yi’ 


where 


( 2 . 16 ) 


( 2 . 17 ) 


( 2 . 18 ) 


( 2 . 19 ) 
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o = 9u] 3v? 3w°3w° dw^aw^ 

^ x>1 8y 3x dx 9y 3x dy 

o _ dw°5w° 

Yxy =" a7'a7 ' 

The curvatures can be handled in like fashion. Specifically: 


which yields 


In the same manner, 


and 


K +AK = 



o o 
= K * +EK ,, 


.0 o 

AK = £K 


K + AK , = 


0 _ d' 


dy 

aV 

3y 2 


{w +ewj) 


-e 


-.2 o 

O Wj 

9y 2 


= K y + £K yi 


and finally 


with 


A Ky = EK yi , 


-.2 

0 . o * a , o o N 

K *y + AK xy = - 2 T^T-( W +£w |) 


dxdy 


3 2 o o 

w n dw, 
= - 2 — 2e- L 


5x5y 5x3y 


O 0 

= K xy + eK *y, 


(2.20) 


( 2 . 21 ) 


( 2 . 22 ) 


(2.23) 


(2.24) 


(2.25) 
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AK *y = £K xy, 


(2.26) 


where 


K =- 


3 2 w° 
3x 2 


->2 o 

° _ dw, (2.27) 

y. ^ 2 

3y 

d 2 o 

„ _ . 'Jh. . 

xy > 9x3y 

The constitutive equation relating the inplane force and out-of-plane bending moment resultants to the 
inplane reference surface strains and the out-of-plane reference surface curvatures is well known and is 
given by 






■a„ 

A 12 

A, 6 

B„ 

B,2 

B,6 _ 

N y 


A 12 

a 22 

A26 

B 12 

®22 

B26 

N *y 


A 16 

A26 

A66 

B 16 

«26 

»66 



B„ 

B 12 

B,6 

D U 

D 12 

D ,6 

M v 


B 12 

b 22 

B 26 

D.2 

d 22 

°26 

y 


B.6 

B 26 

B 6 6 

D,6 

D 2 6 

&«L 


o 

K 




(2.28) 


The Ay, By, and Dy are defined as follows: 

H/2 

J 


H/2 


H/2 


A u = 


Qijdz 


B .j 


1 


QijZdz 


D .j = 


-H/2 


-H/2 


J 


-H/2 


QijZ 2 dz ■ 


(2.29) 


The Qjj are the reduced ply stiffnesses, defined for a state of plane stress, transformed into global 
coordinates. The A matrix relates the force resultants to the reference surface strains; the D matrix relates the 
bending moment resultants to the reference surface curvatures; while the B matrix couples the force 
resultants with the reference curvatures and couples the bending moment resultants with the reference 
surface strains. The variations in the stress resultants in terms of the variations in strains and curvatures 
become 
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AN x ' 


A. 

a ]2 

A .6 

B„ 

B,2 

B 16~ 


AE x I 

AN y 


A,2 

a 22 

A 26 

8 12 

®22 

B 26 


A S 

AN xy 


A 16 

A 26 

A 66 

B ,6 

B 2 6 

B 6 6 


A ?xy 

AM x 


B„ 

B,2 

B 16 

D„ 

D, 2 

D,6 


A < 

AM 


B 12 

B 22 

®26 

D.2 

D 22 

B>26 


AK y 

y 

AM > 


_b 16 

»26 

B 66 

D,6 

D 26 

C>66. 


Ax* 0 ' 


For symmetric laminates the B matrix is zero, providing no coupling between the inplane and out-of-plane 
responses of the plate. Assuming symmetric laminates and substituting the expressions from Eqs. (2.13), 
(2.16), (2.19) into Eq. (2.30) yields the following expression for AN x : 


AN,- A „(ee* +e'e^) + A l2 (ee° y + ee° y ) +A l6 (ef, y +ty° xy ) 

= e(A,,e x + A 12 e y + A ]6 y xy ) +e (^,,6^ + A 12 e y ^ + A 16 Y xy2 ) 


(2.31) 


= eN + £ N > 
*1 *2 


where 


r° 


N x, * A l\\ + A 12 £ yi + A 16% 

N x 2 = A 11 £ x 2 + A 12 £ y 2 + A 16Y X y 2 . 


(2.32) 


Similar expressions are obtained in Eqs. (2.33) through (2.42) for AN , AN . AM , AM , and AM , 

y *y * y xy 

namely, 


AN v = A 12 (ee“ + e 2 e x ) + A„ (ee° + eV ) + A 26 (et^ + eV v ) 


2 o 


£ ( A 12 £ x, + A 22 £ y, + A 26^xy,) +£ ( A 12 £ x, + A 22 £ y 2 + A 26Y X y„) 


* x y 2 


(2.33) 


AN y = eN yi + eN y2 


AN xy = A 16 (be; + b 2 b: 2 ) + A 26 <££* + e 2 e;) + A 66 (B Y : yi + bV x>2 ) 

= e ( A 16 b x ^ + A 26 e y[ + A 66 y xy ^) + e (A 16 e^ + A 26 e y ^ + A^y^) 


^ y = £N yi l£N y , 


(2.34) 
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AM x = D 11 e< i + D 12 eK;+D 16 £< yi 
= e(D n K° + D 12 k“ + D, 6 k“ ) 


AM. = eM. , 


(2.35) 


AM y = D, 2 ek^ + D 22 ek^ + D^^xy, 

= £ (Dj 2 K x ^ + D^Ky^ + 

AM y = £M y/ 


(2.36) 


AM xy = D 16 EK x, + D 26 EK y, + D 66 EK xy, 


= E (D 16 K; + D 26 K; + D 66 < y ) 


AM xy = eM xy; 


where 


Ny, - A 12 E x, + A 22 E y, + ^■26^iy l 

N y2 = A 12 e; + A 22 £; j + A 26 Y ° y2 

Nxy^A.^ + A^ + A^ 

Nxy 2 = A 16 e, 2 + A 26 e y2 + \(fxy 2 


K ) = ^n< i + D n< / , + D A l 


(2.37) 


(2.38) 


(2.39) 


(2.40) 


My, = D 12 k 


°22 y. 


D 26 <y, 


(2.41) 


M xy, - D 16 K x. +D 26 K y, 


^ 66 K x y ,- 


(2.42) 


As stated previously, interest centers on the variation of the total potential energy with respect to variations 
in the displacement field. That is, 
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An = (n + An) -n. 


(2.43) 


Substituting Eqs. (2.6) and (2.7) into Eq. (2.43) yields the desired result, which can be written as 


An = \ If { [ N * E * +N * Ae " + ' AN x< +AN ^3 


+ [Vy + N yAe° + AN y eJ + AN y A £ ;] 

+ [n^+n^+an^+an, XJ 

+ [ M x K x + M x AlC x + AM x K x + AM x A k °J (2.44) 

+ [M yK ; + M y AK y + AM y K° + AM y AK y ] 

+ [m„<, + M„A<, + AM„«; y + AM„A<J 
- [ N .< + N ,<, + N,X, + M,< + M (K ; + M 1( <J } dA 

+ ^load ( U ° + Au °- V ° + Av °- w ° + Aw°) - n |oad ( U°, V°, W°) - 

where, until more is said about the specific loading, n, oad must remain in symbolic form. Expanding Eqs. 
(2.44) according to Eqs. (2.12), (2.15), (2.18), (2.21), (2.23), (2.25), (2.31), (2.33), (2.34), (2.35), (2.36), and 
(2.37) yields 

An = { N x(< + e 2 <) + ( £ N ) , i + e 2 N ]lj )e“+ (eN x + e 2 N X2 ) «, + £2 <) 


+ N y (£ E ° +e 2 e!) + (eN + e 2 N ) e ° x + (eN +e 2 N y ) (ee° + e 2 e“) 


+ N xy Ky, + £2 Yxy 2 ) + (eN,^ + E^) Y ° y 

+ (eN xy +e 2 N xyi ) (ey xyi + e 2 Y ° yj ) 


(2.45) 


+ Mek° + eM ¥ k + eM r £k° 

” X Xj X | 


+ M £K° + EM K° + EM £K° 
y j\ j\ y j\ j] 

+ M xy £K xy, + EM xy,<y + ^xy^xy, > <« A 
+ n ioad ( U° + eu°, V° + £V°, W° + EW°) - n, oad ( U°, V°, W°) ' 


Combining terms with like powers of e produces 
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AH = \ J J { £ [ N x e x, + + N yE ; + n Yi e; + N iyY : yi + N xy) Y° y 

A 

+ M x K x, + M X,< + M y K y, + M y, K y + M xy K xy, + M xy, K xy 1 
+ [ N x < 2 + + N x, e x, + N y E y 2 + N y 2 E y + N y , E y, 

+ + + ] 
+ [ N ., £ ", + N .,< + VJ, + Vn + N -» A , + N -y,< J 

+ e 4 [N x x;H.N yi e; i+ N <A j )dA 
+ n to.d <“° + «•!. v” + e<, w° + EW°) -n loia (»”. .”,*•) 
Equation (2.46) can be expressed as 


( 2 . 46 ) 


An = eOj + e n 2 + e n 3 + e n 4 , (2.47) 

where is the first variation, which is linear with respect to the variations in the displacements, E n 2 is 
the second variation, which is quadratic with respect to the variations in the displacements, etc. 

If the displacement field (u°, v°, w°) is a solution that renders the total potential energy n of a structure to 
be at a minimum, then the change in total potential energy An must be positive for all kinematically 
admissible variations in the displacement field (£u°, £v°, £w°) . The first variation of total potential energy 
is sign-dependent because it is linear with respect to displacement variations. Hence it is possible to have a 
decrease in total potential energy for at least one variation of the displacement field. Therefore, a necessary 
condition for the total potential energy of a structure to have a minimum value in the neighborhood of an 
equilibrium displacement field is for the first variation, ell, , to be zero for all kinematically admissible 
variations in the displacement field. Setting £llj equal to zero leads to the equilibrium conditions for the 
structure. 


The next higher order term, £ 2 n 2 , must be investigated to determine if a minimum of the total potential 
energy exists (as opposed to a maximum or an inflection point). The second variation of total potential 
energy is quadratic with respect to the variations in the displacements. Thus it is sufficient to use this term in 
determining the conditions for which a minimum exists. If £ n 2 equals zero for even one nontrivial 
variation in the displacement field, then £ 3 n 3 must be investigated. 
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The second variation of total potential energy of the plate is a function of the displacement field of the plate. 
The displacement field is, in turn, a function of a loading parameter p. When p = 0, the plate is in stable 
equilibrium in its unloaded state and the equilibrium displacement field has the trivial solution 
u° = v° = w° = 0. Increasing p by sufficiently small amounts will normally yield nontrivial equilibrium 
solutions such that £ Il 2 > 0 (positive definite). The loading parameter p is said to have a critical value p cr 
at the smallest value of p that renders e 2 I1 2 > 0 (positive semi-definite). In general, determining the 
conditions which cause a functional to change from positive definite to positive semi -definite can be 
accomplished using calculus of variations by setting the first variation of the functional equal to zero. 
Accordingly, the value of p cr can be determined by setting the first variation of the second variation of total 
potential energy equal to zero. This approach is known as the Trefftz stability criterion. 

Expanding the second variation of the stress resultants in the second variation of the total potential energy by 
using Eqs. (2.32), (2.38), and (2.39) produces 


^2 2JJ (^11 £ x 2 + ^12 e y 2 + ^l6Yxy 2 ) 

A 

+ (Ny £ y 2 + (A 12 £ x 2 + A 22 £ y 2 + A 26 Y° y ^)£y) 

+ { N xyYxy 2 + ^ A 16 £ x 2 + A 26 £ y 2 + A 66^xy 2 ) X 


.0 \ * r o 
*y 




Using the definitions of Eq. (2.28), (2.48) becomes 


(2.48) 


! n 2 = i|J ( 2N,e; + 2N y e; 2 + 2N„^ yi + N I £^ + N,^ + N,^, 


+ M ..< + M , I < I + M .,,<,, HA 


(2.49) 


The variations in strains and curvature are now expanded by Eqs. (2.14), (2.17), (2.20), and (2.27) to 
produce 
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(2.50) 


dy dy dy 

f du? dv? dw° dw° dw° dw? ^ 

xy > ' dy dx dx dy dx dy ' 

+ M x K° + M K° + M K° } dA • 

For initially flat plates w° = 0. For this analysis, it is assumed that the plate is loaded inplane to the 
buckling condition by specifying the inplane force resultants N x , Ny, and N X y Furthermore, to study stability 
with respect to out-of-plane displacements, it is necessary that the variation occur only in the out-of-plane 
component, w°. These conditions imply the following relations are true: 


dw° _ dw° _ o 
dx dy 

o o rt 
u i = v i = ° 


e o = <K = 0 (2.51) 

*> dx 

e°=^=o 

y ' dy 

= ^ + ^= 0 . 
r * y ' dx dy 

where, as indicated in the 3rd relation, the subscript 1 has been dropped from the out-of-plane displacement 
to simplify notation. Using Eq. (2.51), Eq. (2.50) reduces to 




+1 


fdw°V fdw°dw 0N ) 

N A"^tJ + xy l dx dv J 


dx dy 


(2.52) 


+ Mx 1 < 1 + M y , K y, + M ,y 1 <y,> dA 


Expanding the variations in the moments by using Eqs. (2.40) - (2.42) yields 
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! nj =iff{N f— )^N f— f + 2N„f— — 1 

2 2JJ ‘V3x ) y V3y ) xy V3x 3y ) 


+ (D n < + D I2 K; i + D 16 < yi )< i 
+ (Di 2 < i + D 22 K° ) + D 26 < y ,)K“ 

+ ^!6 K x, + ^26 K y, + ^ > 66 K xy | ) K xy | } ' 

Substituting Eq. (2.27) for the variations in the curvatures (no subscripts on w) results in 

2 2 J J X ^3 x J Aa y ) xy v 3 x 3 y ) 


+ D, 


f-,7 o\ 

a w 


Idx 2 


+ 2D 


Vw-lfsV 1 


+ 4D 


16 


'V 

d w 

j 


12 

faVM a n fdVMfaVM ja 

La y J + 26 [a y 2 J Uxdy J } 


dx 2 


{dy 

f^o 


faVY fa 

Mv I +4D X- 


w 


x5y 


Integrating the coefficient of 4D 66 by parts twice using Green’s theorem yields 


(2.53) 


(2.54) 




r a 2 <a 

a w 

Kdx 2 ) 


a 2 o'i 
a w 

Uy 2 j 


8A 


Figure 2.1 depicts an arbitrary plate having area A enclosed by a smooth boundary 3A with an inplane 
outward normal vector n . The outward normal vector is composed of the components n and n as used in 

x y 

the boundary integral of Eq. (2.55). The term ds appearing in Eq. (2.55) is an infinitesimal arc length along 
boundary 5A traversing the plate in a counterclockwise direction. 
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z 



Figure 2.1 Arbitrary Plate Geometry 


For polygonal plates which have w = 0 along the boundary, the boundary integral appearing in Eq. (2.55) 
is zero (see appendix A). Making this substitution into (2.54) yields 



2.5 Nondimensional Form 

The development to this point has been for a polygonal plate of arbitrary geometry. The second variation of 
total potential energy is now nondimensionalized for the case of a plate with a planform geometry of an 
isosceles trapezoid. This geometry was first presented as Fig. 1.1, and is repeated here for convenience in the 
nondimensional coordinate system. 
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Figure 2.2 Trapezoidal Plate Geometry 


Referring to Fig. 2.2, the geometric nondimensional parameters aspect ratio and taper ratio are defined as 


aspect ratio 

A.R. = — 

w 2 



w 2 

(2.57) 

taper ratio 

T.R. = — • 
W, 


The x and y coordinates are nondimensionalized according to the following relationships: 


«-E 

„ y 

(2.58) 


2.5 Nondimensional Form 22 



Chapter 2: Development of Analysis 


The differential area dA appearing in Eq, (2.56) equals dxdy in the coordinate system of Figure 2.2 which, 
when nondimensionalized, equals Lbjd£,dr| . Using the chain rule of differentiation for the derivatives of w 
yields the relationships 


3w° 

1 3w° 

aV 

dx 2 

i aV 
L 2 a^ 2 

5x 

~Ldl 

3 2 0 
d w 

i a V 

3w° 

1 3w° 

ay 2 

~ b, 2 an 2 

ay 

b,ar) 

3 2 0 

0 w 

i aV 


3x9y Lb,3^dr| 


Substituting Eqs. (2.59) into Eq. (2.56) yields 



(2.59) 


(2.60) 


Equation (2.60) can now be divided into two parts, one dependent only on the bending stiffnesses Du, D 12 , 
etc., and the other dependent only on the inplane force resultants. The division takes the form 

E 2 n 2 = u B + u ls , (2.61) 


where 



2 (D]2 + 2D&6) (a 2 w° 


L 2 b 2 


4D, 


(d_ W° \ { ~”~'26 

' Lb 2 



(2.62) 
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and 


U 


is 


1 f f f N * P W °Y N y pW°V 2N xy pw°V 3 w°'| 


(2.63) 


To nondimensionalize the material properties in the fashion of Nemeth [7], Eqs. (2.62) and (2.63) are first 


multiplied by 


Lb 


: , resulting in 


Vd„d 22 

Lt* 1 i r r j p ii 

7d„D 22 B 2Jj{lW D 22 


^2 oV 
d w 


2 (^i 2 + ^ 66 ) d 2 \v° 

U* J ' X D 22 


Vw°l l_ 2 E 

^T ! 2 J + b 2 V D ! 


r.2 0 v 

d w 

.in*; 


b, |D n D I6 

+ 4 r ,l 


1^224/n 3 D 


d w 


o 0 

1 ' 3 W 


*■) 


L ^22 ^26 

"b 4 I 4/ 


K 2 Jwn 


/ a 2 oW -,2 

d w 


|D "<k 


D 


22 


: m >dWn 

aV 


Lb, 


Td^d 


22 


= lff{ N,b 2 

2 J A JU 2 7 ^ b;v D „**D a 

b W D n n ^ DlI D^ ^ ^ ^ 


(2.64) 


The nondimensional coefficients appearing in Eq. (2.64) can be defined: 


b, Pu 

a = —4 — 

MD 22 


P = 


D 12 + ^66 
7^11^22 


D 


Y = 


6 = 


16 




22 


D 


26 




k. = 


N.b 


71 


Nb 


k y = “z 


TlTD. 


22 

N xy b^ 


7t 2 t/D^ 


(2.65) 


This results in nondimensional forms of the two components for the second variation of total potential 
energy. These components are 
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2.6 First Variation of Second Variation of Total Potential Energy 

The nondimensional functional for the second variation of total potential energy is given by 

fl 2 = U B + U Is . (2.68) 

where U B and U IS are given by Eqs. (2.66) and (2.67). The buckling condition is determined by setting the 
variation of this functional equal to zero, i.e., 

sfl 2 = 8 U b + 8U is = a (2-69) 

This equation will be referred to as the stability equation. Interchange of the variational and differential 
operators and standard steps in variational calculus provides the following relationships: 


8 

>| 

8 



Un 2 


3 2 w°9 2 5w° 

2 T ~ 

d$ 2 d? 

^d 2 w°d 2 8w 0 
3t| 2 dn 2 

8 2 w°8 2 8w° afow 0 8 w° 

8^ 2 8£ 2 + dt, 2 8q 2 

a V ^8w° #8w° 8w° 

d$ 2 ^ 8^ 2 

8 2 w°8 2 5w° 8 2 8w° afow° 

an 2 8^8q + dn 2 


(2.70) 
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Applying Eq. (2.70) to Eqs. (2.66) and (2.67) yields 


and 


8U b 


JJ 


2 d 2 w°3 2 8w° 

{a 

d£ 2 d% 2 


+ P 


5 r ^ O 0 

w crow 

3? 


~\2 q o 2 o 

dOW 0 W 

~Ww- 


1 3 2 w°9 2 8w° 

a 2 3t) 2 dr) 2 


+ 2ay 


28 

+ — 


9 2 o -v2j o 

w crow 

L^ 2 

f9 2 w 0 3 2 8w° 


« Lari 2 


-v2o o -\2 o 

crow o w 

a^ 2 d^. 
a 2 8w° aV 
+ an 2 


d£dr) > 


(2.71) 


8U IS = 



+ a IhT 


+ 


58w°3w° 

% an 



(2.72) 


Note that these functionals are symmetric and bilinear with respect to w and 5w. 


2.7 Out-of-Plane Displacement Approximation 

At this point, it is necessary to know w and 8w as functions of £ and r|. The Rayleigh-Ritz method assumes 
w(x,y) to be an infinite series of the form 


w(x,y) = XX W ijM x) flj ( x, y) . (2-73) 

i = lj = I 


The Wjj are to-be-determined weighting functions; (Jr and are functions that must match the kinematic 
(displacement) boundary conditions at the ends and side edges of the plate, respectively. The kinematic 
boundary condition to be met for a simply supported boundary is that the displacement must be equal to zero 
at the boundary. The displacement and the slope, or derivative of the displacement, must equal zero to meet 
the kinematic boundary condition at a clamped boundary. The bending-twisting coupling terms (D 16 and 
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D 26 ) included in this analysis cause the buckled mode shape to be skewed (like a rectangular plate under a 
shear load) if these terms are of the same order of magnitude as D n , D 22 , D 12 , and D 66 . This skewed mode 
shape can be represented by combinations of symmetric and antisymmetric functions of x and y in Eq. 
(2.73). Because of this, the functional form of <j>j and £2j must contain complete sets of symmetric and 
antisymmetric functions of x and y. For a plate having the geometry of Fig. 2.2, the series along the length 
the plate is very straightforward, namely, 


(Jr = sin 
<J>j = cos 


17CX 

_17 


V ^X 

T 


- cos 


(i+1) 


JCXl 

TJ . 


simple supports 
clamped supports. 


= 1 , 2 ,... 


(2.74) 


These functions contain the complete set of symmetric and antisymmetric functions. The functional forms in 
y, across the width of the plate, however, are not as straightforward because the width of the plate varies as a 
function of x. Thus, the £1 , as indicated in Eq. (2.73) must be functions of both x and y. To develop the 
functional form of the (x, y) , consider first a more simple geometry shown in Fig. 2.3. 


y 



Figure 2.3 Simple Geometry 
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In this figure, h (x) represents a variation in the geometry. For simple supports along AB and CD, an 
admissible set of functions is 


sin 


jgy ~ 

-h(x)_ 


(2.75) 


1 j = 1, 2, ... 


This same function can be used for the isosceles trapezoid of Fig. 2.2 by introducing the simple coordinate 
transformation 


y = y + b(x) , 


(2.76) 


where 


b(x) = 

as illustrated in Figure 2.4. 

y 



(2.77) 


Figure 2.4 Coordinate Transformation Illustration 
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Equation(2.75) then becomes 


sin 


jtty “ 

— > sin 

'jjc(y + b(x)) 

_h (x)_ 


L 2b (x) J 


= sin 


jrcy 

l_2b (x) 


Jn" 

2 . 


■jjcl f jTiy “I 

= sin — cos + 

L 2 J L2b(x)J 


cos 


EH 


fry 

2b (x)_ 


(2.78) 


Investigating individual terms of Eq. 2.77 for specific values of j reveals a pattern identifying the sets of 
symmetric and antisymmetric basis functions. Consider the following terms: 


j = 1 — » cos 
j = 2 — » -sin 
j = 3 — > -cos 
j = 4 — > sin 


ny 


L2b(x)J 

ny 

-b(x)_ 

3rcy 

_2bOOJ 
2?cy " 
bOO- 


first symmetric function 


first antisymmetric function 


second symmetric function 


second antisymmetric function. 


The complete set of symmetric and antisymmetric functions is thus, 


symmetric set: 


antisymmetric set: 


1) 


>] 


3S ( 2 j — L . . 

L 2b(x)J| j=1 

. r *y " 

in j 

Lb(x). 


(2.79) 


(2.80) 


i = 1 . 2 .... 


The kinematically admissible form for w(x,y) can now be written as 


w(x,y) = 

i = ij = i 


W^cos 


(2j- 


1) 


icy 

2b(x)J 


+ W°sin 


7ty 


L b (x) 




(2.81) 


where and W° correspond to the to-be-determined constants associated with the even and odd, i.e., 
symmetric and antisymmetric, respectively, functions of x and y. Remember that (Jr (x) is the complete set 
of symmetric and antisymmetric functions of x, and takes the different forms of Eq. (2.74) depending on the 
boundary conditions at the ends of the plate. 
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It is now necessary to write the out-of-plane displacement approximation in terms of the nondimensional 
variables according to Eq. (2.58). Using the definition of b(x), 


where 



_ n , 

1 + el; 


(2.82) 


e = 



(2.83) 


Equation (2.81) can then be written as 


w 


(5.11) = £XM5)[w*cos[(2j-l)? T ^ +< sin 


nt| 

1 1 + e^J 


Making the following substitutions, 


(2.84) 


¥j (fj) = cos (f) 

(gj) = sinfgj) , 


(2.85) 


yields the series approximation for the buckled shape of the plate in terms of the nondimensional variables, 
namely, 


"(M = XX <t, i ( ^ )[W b' ,, j [f j ( ^ Tl)1+W ? r j [ 8i ( ^ 1 l )]] • 

l = lj = l 


( 2 . 86 ) 


The variational of w, 5w, has the same functional form, specifically 
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8w( 5,T1) = X [5W pq V q [f q ( ^ n)1 +5W ?q r q lg q ( ^’ ^ )]] ' 

p = ] q = 1 


where the series indices are chosen to be different so the product of terms in Eqs. (2.71) and (2.72), which 
are now represented by series, are properly computed. 

2.8 Derivatives of Out-of-Plane Displacement 

To proceed further it is necessary to compute the first and second derivatives of w (^, T|) and 8w (£, T|) 
appearing in U B and U IS of Eqs. (2.71) and (2.72). These derivatives are: 


where 


* -XI< c iX +c iX> • 


N N 




(2.88) 


(2.89) 


a 2 o N N 
d w 




= XX{C^ + C»w»} 


i = lj = l 


(2.90) 


where 
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(2.91) 


dw° 

3ti 
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y.y{N>5+NiX). 


i = ij = 1 


(2.92) 
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where 


where 


where 



(2.93) 


(2.94) 



35w° 
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y y {c, ie 
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5W e +c 10 5W° } ; 
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N N 
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(2.95) 


(2-96) 


(2.97) 


(2.98) 


(2.99) 
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N N 


^ p = Iq = 1 


^ = ii (N“SW^N“6W« } ; 
^ p= , q=1 


and 


d^w 0 


N N 




p = 1 q = 1 


( 2 . 100 ) 


( 2 . 101 ) 


(2.102) 


2.9 Matrix Form of the Stability Equations 

Substituting Eqs. (2.88), (2.90), (2.92), (2.94), (2.96), (2.98) - (2.102) for individual terms appearing in the 
equations for U B and U IS yields the following important relationships: 
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(2.105) 
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faw°ua6w 0N j fa8w ON j raw 0N 

N N N N 

k Y X X y { t ( C !, E Ka + N ii EC oa) WE + (C‘°N 1E + N 10 C 1E ) W°] SW E 

S / j / j X a X j 1 1 ' ij pq li dq ' i? ' ii na n nn ' n * nr 


i = 1 j = 1 p = 1 q = 1 


ij pq 


u pq u pq 7 


pq 


+ [ (C.^n'^ + N^C 10 ) W E + (C 10 N 10 + N I0 C 10 ) W°]8W° } 

v U pq ij pq / ij v ij pq ij pq ’ ij J pq 1 ' 


(2.110) 


Substituting Eqs. (2.103) - (2.107) into Eq. (2.71) and assuming that the integrations can be performed on a 
term-by- term basis, 5U B can be expressed as follows: 
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( 2 . 111 ) 
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where 
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(2.115) 


Substituting Eqs. (2.108) - (2.110) into (2.72) and assuming that the integrations can be performed on a 
term-by-term basis, 5U IS can be expressed as follows: 
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pq 


so, s - £i £ i t (k 0 e : ,w e + k - w°) 8w; 
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where 
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K, 


EE 

Jijpq 




(2.117) 


0 -b($) 


K, 


1 WO 

OE = r n 2 v r-lO^lE^ n ,, kt 10 vt 1E , n ,, /r .10 kT IE . xt io^ie 

5 ijpq J J 


>*K** ; 


(2.118) 


0 -b(5) 


K, 


EO 

Jijpq 


-/ } <*v; e c + ^ 


(2.119) 
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(2.120) 


Combining Eqs. (2.111) and (2.116) yields the series form of the stability equation as 

N N N N 


/ E 

pq 


8n = SI 1 1 ([<^ w '-^L w u> + (< w "-iC w “,]6w; 

i = 1 j = 1 p = 1 q = 1 

+ [< k E, w 5 -p k S 1 w & - (■M-pkSLwJ)]*^ = o 


( 2 . 121 ) 


The p appearing in the equation is the loading parameter associated with the inplane force resultants which 
can be factored out of the geometric stiffness matrix. Note that a negative sign is factored out with p. The 
load factor p is defined such that in Eq. (2.65), 

N* = PN X N y = pN y N xy = pN yx , (2.122) 

where N x , N y , and N xy are the inplane force resultants caused by a load much smaller than the buckling 
load. These quantities are referred to as the prebuckling force resultants. Buckling occurs when N x , N y , and 
N xy are multiplied by p cr . Since the variational displacements 6wj^ and 8W° q are arbitrary, their 
coefficients must vanish. This results in an eigenvalue problem for loading p = p cr that can be expressed in 
matrix form as 
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The terms K? e , Kz° , KZ E , k 2° are elements of the so-called total stiffness matrix [K T ] and the terms 

j ljpq 1 MM hjpq 1 »Jpq 

K EE , K E0 , K° E , K°° are elements of the so-called geometric stiffness matrix [Kq] where 

^ijpq ^'jps Mpq 
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K. c I K g 


(2.124) 


See appendix B for further detail in developing [Kp] and [K G ] . 

2.10 Prebuckling Inplane Force Resultant Distribution 

As stated previously, the equilibrium equations are obtained by setting the first variation of total potential 
energy equal to zero. The resulting equilibrium equations are then solved for the spatial distributions of the 
force resultants N x , Ny, and N xy . In order to simplify this analysis, however, a relatively simple prebuckling 
inplane force resultant distribution is assumed for the plate rather than solving the equilibrium equations. 
The prebuckling inplane force resultant distribution is assumed to be caused by a known shortening of the 
plate in the x direction and is of the form 


N 

x W(x) 

N y = 0 


(2.125) 


In Eqs. (2.125), P is the applied load required on the ends of the plate to produce a given end shortening, and 
W(x) is the width of the plate, the width varying linearly with x (see Fig. 2.2). This simplification of the 
force resultant distribution is based on finite element results obtained using the general purpose finite 
element program ABAQUS [9] showing that N x is significantly greater than N y and N xy for the range of 
geometries under consideration. ABAQUS is well suited to analysis of composite structures because it 
handles anisotropic material properties quite easily. In this case, the ABD matrix of a particular stacking 
sequence is computed using a separate analysis and is entered as input data to ABAQUS. Table 2. 1 presents 
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maximum values for N y and N xy , normalized by the maximum value of N x , for various geometries and 
stacking sequences.The results in Table 2.1 are taken from the finite element studies. The material properties 
used here and throughout are given in Table 2.2. The graphite/epoxy properties represent AS4/3502 [7], the 
material used in the experiments. Aluminum is used to represent the isotropic case and these properties are 
also given in Table 2.2. 


Table 2.1 Maximum Inplane Force Resultants for Various Cases 



1.0 1.0 Isotropic 0 0 

1.0 1.0 [+30 6 /-30 6 ] s 0 0 


1.0 1.0 Isotropic 0 0 


1.0 1.0 [+30 6 /-30 6 ] s 0 0 


1.0 L5 Isotropic 0.084 0,118 


IS) IS 1+ 30«/-30<& 0.075 0.117 


1.0 3.0 Isotropic 0.148 0.213 


1.0 3.0 [+306/-30 6 ] s 0.127 0.216 


s 


2.0 3.0 Isotropic 0.082 0.134 


Note that for a given plate geometry there is very little difference in the ratios of Ny/N x and N xy /N x between 
an isotropic and a highly anisotropic material, indicating that the geometry plays a much stronger role than 
material properties in the force resultant distributions. Also note that for the most severe taper ratio, N x and 
N xy are but 20% the magnitude of N x . 


Table 2.2 Material Properties for AS4/3502 Graphite Epoxy and Aluminum 
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The simple assumption in Eq. (2.125) for N x was made because it agrees fairly well with finite element 
results and is straightforward to incorporate into the computer code written for this analysis. The assumed 
form of N x must be nondimensionalized as well. From Eq. (2.125), 


where 


P 

WOO 


p 



1 + eE, 


(2.126) 



(2.127) 


Figure 2.5 is a plot of the assumed distribution of N x along the length of the plate for various geometries. 
Note that the relation N x vs. x is independent of stacking sequence. 
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Implicit in Eqs. (2.125) is the assumption that N x is uniform across the width of the plate. Figures 2.6 and 
2.7 compare the assumed values for N x with finite element results across the width of the plate at several 
locations along the length for isotropic and [90] 24 composite plates having aspect ratio of 1 .0 and taper ratios 
1 .5 and 3.0, respectively. The solid lines show normalized finite element results, while the dashed lines show 
normalized assumed values for N x , both of which are normalized with respect to N x . The highly orthotropic 
[90] 2 4 laminate is chosen as a contrast to the isotropic plate. Figures 2.8 -2.11 show finite element results, 
again normalized with respect to N x , revealing the spatial distribution of N y and N xy in the plate. 
Remember that these quantities are assumed to equal zero in this analysis. For comparison, Figs. 2.12-2.17 
show results for the same geometries for isotropic and [0] 24 composite plates, plates with an extreme in 
orthotropy. 
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Figure 2.6 Assumed vs. Finite Element Values for : A.R. = 1.0, T.R. - 1.5; Isotropic, [90]24 
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x/L=0.85 



x/L^O.65 


x/L=0.45 


Assumed N xy =0 • 


x/L=0.25 


ABAQUS, Isotropic 
ABAQUS,[90] 24 


Figure 2.10 Finite Element Results for Njy/N* : A.R.=1.0, T.R.-1.5; Isotropic, [90]24 


x/L=0.85 


x/L=0.65 


x/L=0.45 



x/L=0.25 


Assumed N xy =0 


x/L=0.03 


ABAQUS, Isotropic 
ABAQUS, [90] 24 
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Figure 2.11 Finite Element Results for N xy /N x : A.R.=1.0, T.R.=3.0; Isotropic, [90^24 
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Figure 2.14 Finite Element Results for Ny/N* : A.R.=1.0, T.R.=1.5; Isotropic, [O]^ 
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By examining these figures, it is seen that for the isotropic plate N x is least uniform near the narrow end. The 
largest error in assuming N y is zero is at the narrow end. Assuming N xy is zero in an isotropic plate also 
results in the most error at the narrow end. For the [90]24 plate, the error in N x is similar to the isotropic case. 
The eiTor in N y is greater than for the isotropic case and the spatial distribution is different. For the isotropic 
case, the errors in N y are both positive and negative, while for the [90]24 laminate, the errors are all positive. 
For the [90] 24 laminate, the errors in N xy are very similar to the errors in N xy for the isotropic case. For the 
[0] 24 case, the errors in N x are exaggerated relative to the errors in N x for the isotropic case. Oddly enough, 
the errors in N y for the [0] 24 case are considerably less than those of the isotropic case, as is the case for N xy . 
The accuracy, or the lack of it, in the prebuckling force resultant distributions has an impact on the buckling 
predictions. Unfortunately, it is not possible to make general statements regarding whether the error is 
always conservative or always nonconservative. Falsely adding inplane compressive forces to the problem 
will result in the buckling loads being predicted lower than is actually the case. Falsely adding inplane 
tensile forces to the problem will result in the buckling loads being predicted higher than is actually the case. 
In addition, the sign of the falsely assumed shear is important. As can be seen from the figures, the sign of 
the error in representing the inplane force resultant distribution varies from location to location in the panel, 
and it varies with panel orthotropy and, to some degree, panel geometry. Figures 2.6 - 2.17 are considered 
quite important to this study, and will help put the buckling predictions obtained from the semi-analytic 
analysis into context when compared to results obtained from the finite element analysis. 


With the prebuckling stress resultants specified, k x , k y and k xy of Eq. (2.65) can be rewritten as 


k = 


nV 

* l 


n 2 J D,,D 22 (l+e£) 


(2.128) 


k y = 0 
k *y = °- 


These specific expressions greatly simplify the components of the geometric stiffness matrix [K G ]. 
Referring to Eqs. (2.117)-(2.120), 


K, 


EE _ 
Jijpq 


- P 

0 -b(£) 


- 5 — d?c' E }di\dt , 

l+e$ ,J pq 
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(2.132) 


With the geometric stiffness matrix posed in this form, the loading parameter p will then have the form 

N 'b? 

(2.133) 

Note that the negative sign has been removed. Recall, it was factored out in posing the problem in the form 
of Eq. (2.121) 

2.11 Computer Implementation 

The stability equation is posed as a generalized eigenvalue problem having the form 

[K T ] {x} = p[K c ] {x}. (2.134) 

E O 

The eigenvector {x} is an arbitrary arrangement of the unknown constants and , while the loading 
parameter p constitutes the eigenvalue. Obviously, there are multiple eigenvalues and associated 
eigenvectors. The matrices [K T ] and [K^] are symmetric because the functionals 8U B and 8U IS are 
symmetric with respect to w and 8w. The symmetry of the stiffness matrices results in real-valued 
eigenvalues. An additional important property to note is that [K T ] is positive definite because the strain 
energy of a plate is a positive quantity. These special properties of the stiffness matrices allow the use of 
simple eigenvalue extraction routines on the computer, resulting in much shorter run times than general 
eigenvalue extraction routines. 

The elements in [K T ] and [K^] are functions of plate geometry, material properties, and the indices i, j, p, 
and q in the series for w(x,y) and 8w(x,y) (see eqs. 2.85 and 2.86). An important factor in developing the 



2.1 J Computer Implementation 


48 



Chapter 2: Development of Analysis 


computer code is defining the order in which the unknown constants are included in {x} because the i and j 
indices define a specific position within the eigenvector, and hence, a specific column in the stiffness 
matrices. In this analysis, the following pattern is used: 


{x} = {wf 1( wf 2 ,. 

w E w E 


w° w° 

vv lN , w 21 , 


w E 
” 22’ 

W E 

** 2N’ 

W° 
"22* ’* 

w° 

*’ ** 2N’ 


w; 

w: 


Nl» 

O 


Nl» 


W E 

VV N2’ •' 

W E 

w° 

VV N2» •* 

W° 

■’ W NN 


(2.135) 


Once the order of the Wy in the eigenvectors is determined, the structure of the computer program is 
relatively straightforward. Both matrices [K T ] and [K^] are composed of four submatrices each having size 
N 2 x N 2 , namely, ( [K EE ] , [K 0E ] , [K EO ] , [K°°] ) which are be computed concurrently. The 
bookkeeping to keep track of position in the eigenvector is handled by two arrays having length N which 
are constructed as follows: 


{P} = { 1, 1, 1, 2, 2, 2, ..., N, N, ..., N} (2136) 

{Q} = {1,2, 3, ..., N, 1, 2, 3, ...N, ..., 1, 2, 3, N}. 

An outer loop in initiated to step through the N 2 rows of the submatrices. For a specific row, the fixed values 
of p and q are determined from the {P} and {Q} arrays (p=P(k) and q=Q(k) for the kth row) while i and j are 
free indices, each ranging from 1 to N. Once p and q are determined for the specific row, two nested loops 
are initiated to step from 1 to N for both i and j indices, each combination of i and j defining a specific 
column within the row under consideration. After the two nested loops are complete and all of the columns 
are accounted for, the outer loop updates to the next row and the process is then repeated until all rows are 
complete. When the outer loop is finished, the four submatrices comprising [K T ] and [K^] are combined and 
a standard eigenvalue extraction routine is used to solve for the eigenvalues and eigenvectors. The critical 
value of the loading parameter, i.e., the buckling load, is then given by the lowest eigenvalue. See Appendix 
C for a more detailed discussion of the procedure used to construct the submatrices. 

2.12 Convergence and Efficiency of Analysis 

Because a series approximation is used for the out-of-plane displacement of the plate, the accuracy of this 
analysis is dependent on the number of terms used in the approximating series. Table 2.3 presents results 
showing the convergence of this analysis over a variety of geometries for both isotropic and highly 
anisotropic plates with clamped end conditions. In each row of the table analytical results for buckling 
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coefficient are normalized by finite element results for buckling coefficient for a specific material and plate 
geometry. The general purpose finite element code AB AQUS is used. 


Table 2.3 Convergence of Analysis for Clamped End Conditions 
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Note that the solution converges more quickly for the isotropic material than for the composite material 
because more terms are required to accurately model the skewed buckled mode shapes associated with 
highly anisotropic material. For the highly anisotropic case, there is less than 1% difference for successive 
solutions between 5 and 6 terms in the series and the series is assumed to have converged. 

While increasing the number of terms used in the approximating series has a positive effect on the accuracy 
of the solution, it has a profoundly negative effect on the computer run time. For N terms in the 
approximating series, the resulting stiffness matrices have size 2N 2 . Doubling the number of terms in the 
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series has the result of quadrupling the size of the stiffness matrices. This has a significant effect on both the 
number of individual stiffness terms that must be calculated and the eigenvalue extraction routine. Table 2.4 
shows the resulting stiffness matrix size and run times for various numbers of terms in the approximating 
series. The run times are reported for a CRAY YMP. . 

Table 2.4 Matrix Size and Run Time for Number of Terms in Series Approximation 

for Clamped End Conditions 
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2.13 Prebuckling Plate Stiffness 

An additional item of interest in the analysis of structures is the slope of the load vs. end-shortening curve 
prior to the onset of buckling. This quantity is referred to as the prebuckling stiffness of the structure. 
Because the prebuckling inplane force resultant distribution is a function of location along the length of the 
plate for tapered plates (rather than a constant as it is for rectangular plates), the calculations for the 
prebuckling stiffness are more involved than those for rectangular plates. Using the assumed inplane force 
resultant distributions shown in Eq. (2.125) yields the following constitutive equations: 

= A n e x + A 12 £ y 

0 = A 12 e x + A 22 e y (2.137 

0 = A 667xy • 


Solving for N x yields 
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N = 
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(2.138) 


Assuming there is no out-of-plane displacement, substituting in for N x and e* yields 
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Multiplying through by dx yields 
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(2.140) 


which can be integrated along the length of the plate to obtain the following equation for end-shortening as a 
function of load P 
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Note that Eq. (2.141) is undefined for rectangular plates, i.e., W 2 = Wj. For rectangular plates, Eq. (2.141) 
simplifies to 
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(2.142) 


Using Eqs. (2.141) and (2.142) yield the following definitions for the prebuckling stiffness of tapered and 
rectangular plates: 
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(W 2 -W,) 



w, 
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w, = w„ 


(2.143) 


The next chapter presents numerically predicted buckling coefficients and mode shapes for plates with 
simply supported ends. The results based on calculations from the series solution are verified by further 
comparisons with finite element results. The chapter following the next chapter presents similar calculations 
for plates with clamped end conditions. Extremes in laminate arrangement, as well as the case of an isotropic 
plate, are investigated to provide insight into the buckling response of plates with trapezoidal planform. 
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Results for Simply Supported Ends 


3,1 Overview of Simply Supported Boundary Condition 

As mentioned in the previous chapter, the nondimensional out-of-plane displacement of the plate is assumed 
to have the form 


w(5.n) = l w ijVjlfj ( 5. +w?r j [ gj (S, ti)]] f (3-D 

i=lj = ! 

where ty. (£) I. _ x 2 n conta i ns the complete set of symmetric and antisymmetric functions of x which 
must meet the kinematic boundary conditions at the ends of the plate. For simply supported boundary 
conditions, the displacement must equal zero at the boundary. The form used in this analysis to satisfy this 
condition is 


♦,(§) = ( 1 + e£) 2 sin (in£) |“ = , 2 . (3.2) 

The ( 1 + e^) 2 term appearing in Eq. (3.2) is included to factor out identical terms that appear in the 
denominator of the integrand resulting from the integrations with respect to r\ in the expression for the total 
potential energy (see Appendix B). If this term is not removed from the denominator, the resulting 
integrations with respect to J; can not be calculated in closed form. 

3.2 Description of Finite Element Model 

To evaluate the accuracy of the Rayleigh-Ritz analysis, the buckling response of trapezoidal plates is also 
investigated using the general purpose finite element program ABAQUS for a limited number of conditions. 
The finite element model consists of a 10 by 10 mesh of 8-node shell elements having 5 degrees of freedom 
per node. A convergence study shows that a model of this size has converged for the aspect and taper ratios 
considered in this study. Because of the tapered geometry and skewed buckled mode shapes, no symmetry of 
the model can be assumed. 
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3.3 Effect of Fiber Orientation on Nondimensional Parameters 

The material properties that directly affect the stability of laminated plates are the elements of the D matrix. 
The elements of the D matrix are governed by the stacking sequence of the laminate and are accounted for in 
the analysis via the non-dimensional parameters a, p, y, and 5 (see Eq. 2.64). Because a is a function of both 
plate geometry and material properties, to simplify future discussions it is redefined as 

b. 

a = u - 

where, 

a = 


D, 


D 


22 


(3.4) 


Figure 3.1 shows how the four parameters a, P, y, and 8 vary with ply orientation and stacking sequence for 
a typical graphite/epoxy composite. Two stacking sequences, [+0(/"®6Js and [iOlfo, are considered, with 0 
ranging from 0 to 90 degrees. 
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From the definitions in Ch. 2, it is seen that the parameters a and P are functions only of D^, D 12 , D 22 , and 
D 66 . Thus they have the same value for both stacking sequences. The parameters y and 5, on the other hand, 
are functions of D !6 and D 26 . These two terms are greatly different for the two stacking sequences. Hence, 
for certain laminates the bending-twisting coupling terms D 16 and D 26 can have a dramatic effect on the 
buckling load of plates. Typically, plates having plies with a specific orientation grouped together, the [+65/- 
0 6 ] s laminates, have higher D 16 and D 26 terms than plates having plies of the same orientation interspersed, 
the [±0] 6s laminate. These effects are the motivation for studying two different stacking sequences. 


3.4 Comparison of Analysis with Finite Element Results 

This section compares analytical results with finite element results.The first case considered for comparison 
is an isotropic material. Analyzing an isotropic material removes the effects of orthotropy and anisotropy, 
revealing how well the analysis handles the tapered geometry. For the isotropic case the parameters y and 5 
equal zero, |3 equals 1 , while a and £ vary as functions of the aspect and taper ratio. Considering an isotropic 
case, additionally, provides a means of comparison with classical solutions for rectangular plates. In this 
analysis, results are presented in the form of a buckling coefficient, defined as 


k 


c 


pw 2 

n 7 ^ 11^22 


(3.5) 


where P is the applied load, W 2 is the plate width at the wide end as shown in Fig. 2.2, and D n and D 22 are 
plate bending stiffnesses as defined by Eq. (2.29). Note that the buckling coefficient, as defined by Eq. (3.5) 
is slightly different than the loading parameter p cr that appears naturally as a result of nondimensionalizing 
the first variation of the second variation of total potential energy. The results are presented in terms of k^ 
rather than in terms of p cr , because the form of the buckling coefficient is similar to that used by Brush and 
Almroth [10], Timoshenko [11] and Nemeth [7] [8] and allows direct comparison with their results for 
rectangular plates. 

Figure 3.2 presents buckling coefficients as a function of plate aspect ratio for various plate taper ratios for 
an isotropic plate. Analytical results are represented by lines and finite element results are represented by 
symbols. The analytical results presented are for four terms in the series approximation for the out-of-plane 
displacement. Finite element results are computed for a number of aspect ratios between 1 .0 and 2.0, and for 
aspect ratio of 3.0. Taper ratios of 1.0, 1.5, 2.0, and 3.0 are considered. 
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A.R. (L/W 2 ) 

Figure 3.2 Buckling Coefficients for Isotropic Plates: Simply Supported Ends 

As can be seen in Fig. 3.2, there is less than 1% difference between the analysis and the finite element results 
for all of the aspect and taper ratios considered with the finite element analysis, the analytic results 
predicting slightly greater buckling coefficients than the finite element results. The relationship for the 
rectangular plate (T.R.=1.0) matches the classical solution exactly, a festooned shape with a minimum at 
k c= 4 [ 5 ] [ 10 ], This festooned shape is the result of mode shape changes, specifically a change in the number 
of half waves along the length of the plate. As the taper ratio increases, the critical buckling coefficient 
increases and the mode changes occur at smaller aspect ratios. This can be interpreted as the taper creating a 
narrower average plate width, causing a stiffening effect in the plate, and decreasing the critical length of the 
buckling half-wave. Note also that the relationships are smoother for plates having higher taper ratios, there 
not being as much sharpness to the cusps in the festooned shapes. 

Buckling coefficients, as a function of plate aspect ratio, are shown for various plate taper ratios of graphite/ 
epoxy composite plates having stacking sequences [0 ] 24 and [90] 24 , respectively, in Figs. 3.3 and 3.4. These 
stacking sequences are limiting cases for investigating how well the analysis handles orthotropy. The 
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parameters yand 5 equal zero, however, P now equals 0.405 for AS4/3502 graphite/epoxy. The results in 
Figs. 3.3 and 3.4 for [0] 2 4 and [90] 2 4 laminates, respectively, are for five terms in the series for the out-of- 
plane displacement approximation, and show the same trends as the isotropic case, specifically, an increase 
in buckling coefficient due to an increase in the taper ratio of the plate. However, there are greater 
differences between the analytical and finite element results for these composite cases than the isotropic 
case. Again, the analytic results are greater than the finite element results. For the [0] 2 4 laminate, the analytic 
solution yields nonconservative results 12% higher than the finite element solution for taper ratio 3.0, 10% 
higher for taper ratio 2.0, and 5% higher for taper ratio 1.5. For the [90] 24 laminate the analytic solution 
yields results 6%, 4%, and 2% higher, respectively, than finite element results. This closer agreement 
between the analytic and finite element results for the tapered geometries of the [90] 24 laminate compared to 
the [0] 24 laminate can be attributed to the more accurate approximation of the prebuckling inplane force 
resultant distribution for the “softer” laminate, as was shown in Figs. 2.6 - 2.1 7. For both the [0] 24 and [90] 24 
laminates, the analysis again yields exact agreement between the analytic and finite element solutions for 
rectangular plates 



A.R. (L/W 2 ) 


Figure 3.3 Buckling Coefficients for [0]24 Laminate: Simply Supported Ends 
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A.R. (L/W 2 ) 

Figure 3.4 Buckling Coefficients for [90]24 Laminate: Simply Supported Ends 

Buckling coefficients as a function of plate aspect ratio for various plate taper ratios of graphite/epoxy plates 
having stacking sequences of [±30] 6s and [+3(y-30 6 ] sl respectively, are presented in Figs. 3.5 and 3.6. The 
analytic results are shown for a five term displacement approximation, and include the effects of geometry, 
orthotropy, and anisotropy. Notice that the sets of relationships for the two different laminates have nearly 
identical shapes, but the buckling coefficients for the [+30^/-30g] s laminate are approximately 20% lower 
than those for the [±30] 6s laminate. This is a direct result of the large values of D 16 and D 26 for the [+30 6 /- 
30g] s laminate with the grouped plys relative to the values for the [±30] ()S laminate with the dispersed plys. 
Figures 3.5 and 3.6 show excellent agreement between than analytical and finite element results for taper 
ratios up to 1.5. The analytical results are approximately 5% lower than finite element results for taper ratio 
2.0 and approximately 1 1 % lower than finite element results for taper ratio 3.0. Note that for the higher taper 
ratios the analytical results are less than those obtained using finite element analysis for the 30 degree 
laminates, while the analytical results are greater than the finite element results for the 0 and 90 degree 
laminates. 
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A.R. (L/W 2 ) 

Figure 3.5 Buckling Coefficients for [±30] 6s Laminate: Simply Supported Ends 



0 0.5 1 1.5 2 2.5 3 


A.R. (L/W 2 ) 

Figure 3.6 Buckling Coefficients for [+30 6 /-30 6 ] s Laminate: Simply Supported Ends 
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The four preceding figures provide a general measure of accuracy of the analysis over the range of 
geometries considered for plates with simply supported ends. Based on comparison with finite element 
results, it appears safe to say that, for laminates of practical interest, this analysis can predict buckling 
coefficients accurate to within 1 1 % for taper ratio 3.0, 6% for taper ratio 2.0, 3% for taper ratio 1 .5, and 
0.5% for taper ratio 1.0. These errors are significantly smaller at higher aspect ratios and are accurate to 
within 5%, 2%, and 1% for taper ratios 3.0, 2.0, and 1.5 at aspect ratios greater than 1.5. Nothing, in general, 
can be stated regarding whether the errors are conservative or non-conservative. 

Before proceeding, a comment is in order regarding the semi-analytical predictions of the buckling 
coefficient sometimes being greater than the predictions based on finite elements and sometimes being less. 
Generally the semi-analytical predictions with 5 to 10 terms in the series should be greater than the finite 
element predictions. This is because the number of degrees of freedom in the finite element model would be 
larger and the semi-analytical prediction should approach the finite element prediction from above as the 
number of terms in the series is increased. However, this assumes that the prebuckling equilibrium solution 
is the same in the finite element model as in the semi-analytical model. We have seen that this is not the case. 
It is believed that using the assumed simplified form of the prebuckling stress distributions is the primary 
cause of the discrepancies between the analytical and finite element results shown in the previous figures, 
and the lack of consistency in the semi-analytic solution approaching the finite element calculations from 
above. 

3.5 Dimensional Relationships 

This section presents analytical results of buckling coefficient over a wide range of ply orientations, stacking 
sequences, and plate geometries for the case of simply supported ends. Because of the demonstrated 
accuracy of the analysis, comparisons are not made with finite element results. The following figures show 
the effect of ply orientation, stacking sequence, and plate taper ratio for AS4/3502 graphite/epoxy composite 
plates. 

Figure 3.7 presents buckling coefficients as a function of ply orientation for (±0) laminates for several plate 
taper ratios and the two different stacking sequences [±0]6 S and [+0 6 /-06l s . The solid lines are results for the 
[±0] 6s stacking sequence and the dashed lines are results for the [+06/-06l s stacking sequence. Different 
taper ratios are indicated by various symbols, and the aspect ratio is maintained at 1.0 for all cases. For Fig. 
3.7, and similar figures, calculations are performed in 1 degree increments, while symbols are shown only 
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every 15 degrees. Note the significant change in slope that occurs near 0=60 degrees for the relationships 
having T.R. =1.0. 



0 

Figure 3.7 Buckling Coefficients for (£0) Laminates: Simply Supported Ends; A.R. = 1.0 

Again note the significant effect the large D 16 and D 2 6 terms, resulting from the grouped plies, has on the 
buckling load. This effect is least for 0 near 0 and 90 degrees and greatest for 0 between 40 and 50 degrees. 
The ply orientation angles for which the buckling load is maximum are given in Table 3.1 . 

Table 3.1 Ply Orientation Angle 0 for Maximum Buckling Load: Simply Supported Ends 


■' ■ 

y ", 

ip ® 


1.0 

45 

40 

1.5 

46 

44 

2.0 

44 

42 

3.0 

42 

42 


3.5 Dimensional Relationships 


62 




Chapter 3: Results for Simply Supported Ends 


The effect of an increase in plate taper ratio causing an increase in buckling coefficient appears to be most 
significant as the plate taper ratio is first increased from 1 .0, i.e., there is a larger increase in buckling 
coefficient between taper ratios 1.0 and 1.5 than between 2.0 and 3.0. To demonstrate this more clearly, 
buckling coefficients, as a function of taper ratio, for several different graphite/epoxy laminates are shown in 
Fig. 3.8. 



Figure 3.8 Buckling Coefficient vs. Plate Taper Ratio: Simply Supported Ends 

Note that for all of the cases shown in Fig. 3.8, with the notable exception of the [0]24 laminate, the slope of 
the curve is greater near taper ratio 1 .0 and asymptotically approaches a maximum value near taper ratio 3.0. 
The curve for the [0] 2 4 laminate, however, actually has a negative slope near taper ratio 1 .0 and does not 
appear to be converging to a maximum near taper ratio 3.0. 

3.6 Buckled Mode Shapes 

An item of interest, in addition to the critical buckling load, is the buckled mode shape. The buckled mode 
shape can be calculated for a specific eigenvalue by using the elements of the corresponding eigenvector as 
the weighting factors in the assumed series for the out-of-plane deflection of the plate. Figures 3.9 through 
3.20 show buckled mode shapes calculated for a wide range of plate geometries, ply orientations, and 
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stacking sequences by using the eigenvector corresponding to the lowest eigenvalue. The mode shape 
calculated using the lowest eigenvalue is often referred to as the fundamental or primary mode. Comparing 
Figs. 3.9 and 3.10 shows the skewing of the mode shape resulting from the large D 16 and D 2 6 terms of the 
[-i-30^/-30 6 ] s laminate. Figures 3.11 and 3.12 show mode shapes for the same pair of stacking sequences as 
Figs. 3.9 and 3.10 but for plates having tapered geometries. 
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It is clear that the Dj 6 and D 26 have a significant influence on the buckled shape of the plate. Figures 3.13 
and 3.14 show the same plate geometry as Figs. 3.11 and 3. 12 but for [±60]^ and [-k605/-60^] s laminates, 
respectively. Note that these laminates, because they are softer in the loading direction than the (±30) 
laminates, buckle with two half waves in the loading direction. Figures 3.15 through 3.18 contrast the mode 
shapes of a very “stiff' laminate with a very “soft" laminate. Note again the laminate that is softer in the 
loading direction buckles in multiple half waves. Figures 3.19 and 3.20 show the mode shapes for two 
different stacking sequences, one highly anisotropic, of a highly tapered geometry. 

This survey of buckled shapes provides insight into the character of the buckling of tapered plates, and the 
important role of laminate stacking arrangement 

3.7 Nondimensional Relationships 

Up to this point, buckling coefficients have been presented for plates with specific geometries, ply 
orientations, and stacking sequences. This section, however will present a series of generic relationships 
showing the effect of individual nondimensional parameters on the buckling coefficients. This series of 
relationships, based on nondimensional parameters, covers a wide range of material properties and plate 
geometries in relatively few graphs. Reference points showing specific stacking sequences are included for 
an isotropic material and ply orientation angles 6 (almost 0), 45, and 84 (almost 90) degrees. Refer to Fig. 
3.1 for nondimensional parameter values at other ply orientation angles. 

Figures 3.21 - 3.23 show the buckling coefficient as a function of the parameter & for various combinations 
of P, Y> and 8 for aspect ratio 1 .0 and taper ratios 1 .0, 1 .5, and 2.0 respectively. Recall from the definitions of 
Eqs. (2.64), (3.3), and (3.4), that a and p are measures of bending orthotropy, i.e., a is a measure of the 
bending stiffness in the longitudinal direction relative to the bending stiffness in the transverse direction. 
Unfortunately, there is not a corresponding simple physical meaning for the quantity p because it is a 
combination of all of the orthotropic terms in the bending stiffness matrix. It is a measure of the curvature 
effects (anticlastic curvature and twisting curvature) relative to the average bending stiffness. The quantities 
y and 8 provide a measure of bending anisotropy, i.e., and D 26 , respectively. It is clear that an increase in 
the anisotropic parameters Y and 8 causes a decrease in the buckling coefficient. This was shown in previous 
figures comparing stacking sequences with grouped versus interspersed plies and it is also illustrated in these 
nondimensional relations. It is also clear that an increase in the orthotropic parameter P causes an increase in 
the buckling coefficient. The orthotropic parameter a appears to have the smallest effect on the buckling 
coefficient, though for values of a > 0.9, buckling coefficients increase more or less monotonically at a very 
slow rate. The cusps due to a changes in mode shape are noteworthy. As observed before in the dimensional 
plots, the relations become^moother for plates with higher taper ratios. Note the general increase in 
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buckling loads as the taper ratio increases, i.e., comparing Fig. 3.21 for taper ratio 1.0 as opposed to Fig. 
3.23 for taper ratio 2.0. 


[±4514,: a=l.T>, P=2J, 7=8=0.08 
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Figure 3.21 Effects of a, 7 , and 5 for several values of P on Buckling Coefficient: 
Simply Supported Ends; A.R. = 1.0, T.R. = 1.0 
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Figure 3.22 Effects of a , y, and 5 for several values of P on Buckling Coefficient 
Simply Supported Ends; A.R. = 1.0, T.R. = 1.5 
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Figure 3.23 Effects of a, y> and 5 for several values of P on Buckling Coefficient 
Simply Supported Ends; A.R. = 1.0, T.R. = 2.0 
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To demonstrate more clearly the effect of P, Figs. 3.24 - 3.26 show the buckling coefficient as a function of 
the parameter P for various combinations of a, y, and 8 for the same geometries shown in Figs. 3.21 - 3.23. 
Again it is clear that an increase in the anisotropic parameter y and 8 causes a decrease in the buckling 
coefficient. These figures also demonstrate that the parameter P has a larger effect on the buckling coefficient 
than the parameter a, especially at values of & between 0.5 and 1.0. Note that the relations become spaced 
closer together as P increases. This indicates that the decrease in buckling load due to anisotropy decreases 
for higher values of p. As the taper ratio of the plate increases, the sets of curves generated by the three 
different values for a tend to overlap each other more. This implies that as the taper ratio of the plate 
increases, the effect of the plate geometry on the buckling coefficient is as significant as the ratio of 
longitudinal and transverse bending stiffnesses. These areas of overlap are significant from a design 
perspective because they provide flexibility when choosing geometry and material properties to achieve a 
specific buckling coefficient 



Figure 3.24 Effects of p, y, and 8 for several values of a on Buckling Coefficient: 
Simply Supported Ends: A.R. = 1.0, T.R. = 1.0 
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Figure 3.25 Effects of (3, y, and 5 for several values of a on Buckling Coefficient 
Simply Supported Ends: A.R. = 1.0, T.R. = 1.5 


Figure 3.26 Effects of P, y, and 8 for several values of a on Buckling Coefficient 
Simply Supported Ends: A.R. = 1.0, T.R. = 2.0 
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The next chapter presents results for plates with clamped ends. Buckling coefficients as functions of plate 
geometry and laminate arrangement will again be discussed. These results will be compare with the results 
of simply supported plate to provide further insight into the buckling response of plate with tapered 
planform. 
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Chapter 4 

Results for Clamped Ends 


4.1 Overview of Clamped Boundary Condition 

For convenience, the assumed form of the nondimensional out-of-plane displacement of the plate is repeated 
here. The assumed form is 


W (5, 1 \) = XZ <t, i^ )[W u'l , jl f j ( ^ 1 l )]+W ° r j [ 8j^' T l)l]* (41) 

i= lj= l 

where (£) ^ = \ 2 n conta i ns the complete set of symmetric and antisymmetric functions of x which 
must meet the kinematic boundary conditions at the ends of the plate. For clamped boundary conditions, the 
displacement and the slope in the ^ direction, or derivative of the displacement, with respect to must equal 
zero at the boundary. The form used in this analysis to satisfy these conditions is 

<t>j(£) = (l+e£) 2 {cos[(i-l)7c£] -cos[(i + l)ji£]}|~ =I ^ . (4.2) 

The term ( 1 + e^) 2 appearing in Eq. (4.2) is included to factor out identical terms that appear in the 
denominator of the integrand resulting from the integrations with respect to r\ in the expression for the total 
potential energy (see Appendix B). If this term is not removed from the denominator, the resulting integrals 
over ^ can not be calculated in closed-form. 

4.2 Comparison of Analysis with Finite Element Results 

As was done for plates with simply supported ends, this section compares analytical results with finite 
element results from ABAQUS.The first case considered is an isotropic material. As was stated previously, 
analyzing an isotropic material removes the effects of orthotropy and anisotropy, revealing how well the 
analysis handles the tapered geometry for clamped ends and simply supported side edges. Recall for the 
isotropic case the parameters yand 5 equal zero, (3 equals 1, while a and e vary as functions of the plate 
aspect and plate taper ratios. Buckling coefficients as a function of plate aspect ratio for various plate taper 
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ratios for an isotropic plate are shown in Fig. 4. 1 . Analytical results are represented by lines and the finite 
element results are represented by symbols. The analytical results shown are obtained using five terms in the 
series approximation for the out-of-plane displacement approximation. 



A.R. (L/W 2 ) 


Figure 4.1 Buckling Coefficients for Isotropic Plates: Clamped Ends 


There is less than 2 percent difference between the analysis and the finite element results for plate taper 
ratios less than 2. Differences of approximately 4.5% occur at aspect ratios less than 1.3 for taper ratio 3.0. 
Note that for the isotropic clamped case, the analytic results are slightly less than the finite element results. 
This is opposite to the situation for the isotropic simply supported case, Fig. 3.2. Comparisons with Fig 3.2 
show that the clamped boundary conditions have a significant stiffening effect on the plate. Buckling 
coefficients can be twice as large for a plate with clamped ends as those for a plate having the identical 
geometry but with simple supports on the ends. Similar trends are exhibited in Fig 4. 1 that were noted in Fig. 
3.2, specifically, as the taper ratio increases the budding coefficients increase, the mode changes occur at 
smaller aspect ratios, and the relationships become more smooth, the cusps not being as distinct. The 
relationship for the rectangular plate, (T.R. = 1.0) again matches the classical solution exactly. To be noted in 
the clamped case, for aspect ratios between 0.5 and 1.0, the buckling coefficient relations for the various 
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taper ratios cross over each other, reversing the effect of taper ratio. This was not the case for the simply 
supported case in Fig. 3.2 

Buckling coefficients, as a function of plate aspect ratio, for various plate taper ratios of graphite/epoxy 
composite plates having stacking sequences [0] 2 4 and [90] 2 4, respectively, are shown in Figs. 4.2 and 4.3. 
The parameters y and 5 equal zero, however, p now equals 0.405 for AS4/3502 graphite/epoxy material. The 
ply properties for this material are as given in Table 2.2 



A.R. (L/W 2 ) 

Figure 4.2 Buckling Coefficients for [O]^ Laminate: Clamped Ends 
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Figure 4.3 Buckling Coefficients for [90]24 Laminate: Clamped ends 

These relationships for the [90]24 laminate have the same trends as the isotropic case, specifically, an 
increase in buckling coefficient due to an increase in the taper ratio of the plate. Note, however, that Fig. 4.2 
reveals a decrease in the buckling coefficient due to an increase in the plate taper ratio for [ 0]24 laminates 
having plate aspect ratios less than 1.7. As seen in Fig. 3.3, this phenomenon is not exhibited in [0] 2 4 
laminates with the same geometries but with simple supports on the ends. This is a further exaggeration of 
the aforementioned reversal effect seen in Fig. 4.1 relative to Fig. 3.2. 

Note that for the [0] 24 and [90] 2 4 laminates, the buckling coefficients predicted by the analysis are greater 
than the buckling coefficients predicted by the finite element approach. This is similar to the simply 
supported case for [0] 24 and [90] 2 4 laminates. There is better agreement between the analytical and finite 
element results for the [90] 24 laminate than for the [0] 24 laminate. This is consistent with the observation, in 
Ch. 2, that the assumed N x distribution is closer to the actual prebuckling stress resultant distribution for the 
[90] 24 composite plates than for the [0] 24 composite plates (see Figs. 2-6 -2.11 vs. Figs 2.12 - 2.17). For the 
[90] 24 laminate, the analytic solution yields nonconservative results 3% higher than the finite element 
solution for taper ratio 3.0, 1.5% higher for taper ratio 2.0, and less than 1% higher for taper ratio 1.5. For 
plate aspect ratios less than 1.4, the analytic solution for the [0] 24 laminate yields results 9% higher than the 
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finite element solution for taper ratio 3.0, 7% higher for taper ratio 2.0, and 5% higher for taper ratio 1 .5. The 
analysis yields the same solution as finite elements for the rectangular case of both laminates. 

Buckling coefficients as a function of plate aspect ratio for various plate taper ratios of graphite/epoxy plates 
having stacking sequences of [+30]g s and [+306/-30g] s , respectively, are presented in Figs. 4.4 and 4.5. 
These cases include all of the effects of geometry, orthotropy, and anisotropy in the analysis. Notice that the 
sets of relationships for the two different laminates have nearly identical shapes, but the buckling loads for 
the [+3(y-30 6 ] s laminate are approximately 16% lower than those for the [±30] 6s laminate. The same trend 
is visible in Figs. 3.5 and 3.6 for plates having the same stacking sequences but with simply supported ends. 
Figures 4.4 and 4.5 show excellent agreement between the analytical and finite element results for taper 
ratios up to 1.5. The analytical results are approximately 5% lower than finite element results for taper ratio 
2.0 and approximately 1 1 % lower than finite element results for taper ratio 3.0. Note again the reversal of 
trends with respect to taper ratio as the aspect ratio increases past 1.1 - 1.2 
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A.R. (L/W 2 ) 

Figure 4.5 Buckling Coefficients for [+30^/ - 3 O 5 ] s Laminate: Clamped Ends 

The four preceding figures provide a general measure of the accuracy of the analysis over the range of 
geometries studied for plates having clamped ends. Based on comparison with finite element results, it 
appears safe to say that for all laminates of practical interest, the analysis can predict buckling coefficients 
accurate to within 11% for taper ratio 3.0, 7% for taper ratio 2.0, 4% for taper ratio 1.5, and 0.5% for taper 
ratio 1 .0 for plates having simple supports on the side edges and clamped supports on the ends. Note that this 
is the same level of accuracy observed for the plates having simply supported boundary conditions on all 
four edges. 

4.3 Dimensional Relationships 

This section presents analytical results of buckling coefficient over a wide range of ply orientations, stacking 
sequences, and plate geometries for the case of clamped ends. Comparisons are not made with finite element 
results. The following figures show the effect of ply orientation, stacking sequence and plate taper ratio for 
AS4/3502 graphite/epoxy composite plates. Figure 4.6 presents buckling coefficient as a function of ply 
orientation for ( + 0 ) laminates for several plate taper ratios and the two different stacking sequences [rh0] 6s 
and [-h0^/-0 6 ] s . The solid lines are results for the [±0] 6s stacking sequence and the dashed lines are results for 
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the W-eds stacking sequence. Different taper ratios are indicated by symbols, while the aspect rauo is 


maintained at 1.0 for all cases. 



Figure 4.6 Buckling Coefficients for (±9) Laminates: Clamped Ends; A.R.-1.0 

Comparison with Fig 3.7, which shows the same geometries and stacking sequences for (±0) laminates 
having simple supports on the ends, shows that boundary conditions have very significant effects on the 
buckling response of the plate other than causing an increase in the buckling load. A simple comparison 
between Figs. 3.7 and 4.6 is that for the clamped ends the buckling coefficient more or less decreases 
monotonically with increasing fiber orientation angle, while for the simply supported ends the buckling 
coefficient increases with increasing fiber orientation angle and then decreases. However, perhaps the most 
significant characteristic to notice in Fig. 4.6 is that for each of the two stacking sequences there is a unique 
fiber orientation angle at which increasing the plate taper ratio has no effect on the buckling coefficient of 
the plate. This angle will be referred to as 9 crit . For stacking sequences having fiber orientation angle less 
than 0 crit , the buckling coefficient decreases for an increase in taper ratio. If the fiber orientation angle is 
greater than 0 crit , the buckling coefficient increases for an increase in taper ratio. This phenomenon of a 
unique combination of plate aspect ratio and ply orientation, either side of which the effect of plate taper 
ratio has drastically different effects, was referred to in discussion of Figs. 4.2, 4.4, and 4.5. It should be 
noted that 0 crit is different for various plate aspect ratios and stacking sequences. As an example of this, Fig. 
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4.7 presents the buckling coefficient, as a function of ply orientation angle, for the same geometries and 
stacking sequences as shown in Fig. 4.6, but with a plate aspect ratio of 2.0 rather than 1 .0. The relationships 
shown in Fig. 4.7 look more like the relationships of Fig. 3.7. Note that there is no critical angle in the 
stacking sequence for this plate aspect ratio. 



0 

Figure 4,7 Buckling Coefficients for (±0) Laminates: Clamped Ends; A.R.=2.0 

Figure 4.8 presents results for plate aspect ratio 1.5. Note that although the relationships for the two 
laminates cross each other, a clearly defined critical angle does not exist. 
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Figure 4.8 Buckling Coefficients for (£0) Laminates: Clamped Ends; A.R.=1.5 
4.4 Buckled Mode Shapes 

An item of interest, in addition to the critical buckling load, is the buckled mode shape. The buckled mode 
shape can be calculated for a specific eigenvalue by using the elements of the corresponding eigenvector as 
the weighting factors in the assumed series for the out-of-plane deflection of the plate. Figures 4.9 through 
4.20 show buckled mode shapes calculated using the eigenvector corresponding to the lowest eigenvalue for 
a wide range of plate geometries, ply orientations, and stacking sequences. Comparing Figs. 4.9 and 4.10 
shows the skewing of the mode shape resulting from the large D 16 and D 2 6 terms of the [+306/-30 6 ] s 
laminate. Figures 4.1 1 and 4.12 show mode shapes for the same pair of stacking sequences as Figs. 4.9 and 

4.10 but for plates having tapered geometries. Figures 4.13 and 4. 14 show the same plate geometry as Figs. 

4.11 and 4.12 but for [±60] 6s and [^60^/-60^] s laminates, respectively. Figures 4.15 through 4.18 contrast the 
mode shapes of a very “stiff” laminate with a very “soft” laminate. Figures 4.19 and 4.20 show the mode 
shapes for two different stacking sequences, one highly anisotropic, of a highly tapered geometry. 
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Comparing Figs. 3.9 with Figs. 4.9 shows that the clamped end conditions force greater curvature for the 
clamped plate because the plate mode shape is moved in towards the center of the plate, although the 
longitudinal location of maximum out-of-plane displacement is the same for the two different sets of end 
conditions. Figures 3.9 and 4.9 show the same trend. Additionally, as seen when comparing Figs. 3.10 with 
4.10, the skewing effects due to material anisotropy are amplified by the clamped end conditions. 
Comparing Figs. 3.11 - 3.14 with Figs. 4.11 - 4.14 shows that the clamped end conditions cause the peak of 
the out-of-plane displacement to shift toward the narrow end of the plate. Comparing Figs. 3. 1 8 and 4. 1 8 
shows that the simply supported boundary conditions cause the [90] 24 laminate to buckle into three half 
waves with a relatively large area in the narrow end having very little out-of-plane displacement, while the 
clamped boundary conditions cause this laminate the buckle into 4 half waves. 

4.5 Nondimensional Relationships 

As in Ch. 3, this section presents a series of generic relationships showing the effect of individual 
nondimensional parameters on the buckling coefficients. This series of relationships, based on 
nondimensional parameters, covers a wide range of material properties and plate geometries in relatively 
few graphs. When possible, points showing a specific ply orientation or stacking sequence are included for 
reference. Figures 4.21 - 4.23 show the buckling coefficient as a function of the parameter a for various 
combinations of P, y, and 6 for aspect ratio 1.0 and taper ratios 1.0, 1.5, and 2.0 respectively. It is clear that 
an increase in the anisotropic parameters y and 8 causes a decrease in the buckling coefficient. This was also 
shown in previous figures comparing stacking sequences with grouped versus interspersed plies. It is also 
clear that an increase in the orthotropic parameter p causes a decrease in the buckling coefficient. Unlike 
Figs. 3.21 - 3.23 for the case of simply supported ends in Ch. 3, however, the orthotropic parameter a 
appears to have a significant effect on the buckling coefficient. Specifically, the buckling coefficient 
increases almost monotonically with an increase in a. There is one cusp due to a changes in mode shape 
and, as observed before regarding the plots of dimensional relations, the nondimensional relations become 
smoother for plates with higher taper ratios. 
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Figure 4.23 Effects of a, y, and 5 for several values of P on Buckling Coefficient: 

Clamped Ends; A.R. = 1.0, T.R. = 2.0 

To demonstrate the effect of (3, Figs. 4.24 - 4.26 show the buckling coefficient as a function of the parameter 
P for various combinations of a, y, and 8 for the same geometries shown in Figs. 4.21 - 4.23. Again it is 
clear that an increase in the anisotropic parameters y and 8 causes a decrease in the buckling coefficient, 
however, in Figs. 4.24 - 4.26 changes in y or 8 cause significant changes in the curvature of the relations in 
addition to shifting them downward, whereas the curvature does not change significantly due to increased y 
and 8 in Figs. 3.24 - 3.26. As shown in Figs. 4.24 - 4.26, the relations become spaced closer together as P 
increases, indicating the detrimental effects of anisotropy on the buckling coefficient decrease for higher 
values of p. Note that there are significant areas of overlap between the three sets of relations shown in each 
figure, and this area of overlap becomes larger for increased taper ratios. This overlap area implies that a 
designer should have a moderate amount of flexibility when choosing ply orientation on stacking sequence 
to achieve a specific buckling coefficient. 
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Figure 4.26 Effects of p, y, and 5 for several values of a on Buckling Coefficient: 

Clamped Ends; A.R. = 1.0, T.R. = 2.0 

It is clear by comparing Figs, 4.21 - 4.26 with Figs 3.21 - 3.26 that clamped boundary conditions on the ends 
can have a significant effect on the buckling load through interaction with material properties. Specifically, 
the ratio of longitudinal bending stiffness to transverse bending stiffness represented by the parameter a has 
a much more significant effect for clamped plates than it does for simply supported plates. 

With the characteristics of the buckling response of plates with trapezoidal planform having been thoroughly 
discussed for plates having either clamped or simply supported ends, and with insight into the influence of 
the various parameters on the response established, attention turns to the experimental investigation. The 
next chapter describes the experiments conducted as part of this study 
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Chapter 5 

Experimental Approach 


5.1 Overview of Experimental Program 

In order for analytical results to be credible and useful, they should be verified experimentally. To 
complement the numerical predictions, experiments were conducted to measure the buckling characteristics 
of laminated composite plates with a trapezoidal planform area over a range of stacking sequences and plate 
geometries. This chapter presents the plate geometries and stacking sequences investigated; describes the 
specialized test fixtures used to provide the required boundary conditions for the specimens; discusses the 
various instrumentation techniques used to measure end-shortening, inplane strain, out-of-plane deflection, 
and buckled mode shape; and defines the specific experimental procedure used to conduct the experiments 
and reduce the data. 

5.2 Test Specimen Geometries and Stacking Sequences 

The variables of interest for this problem can be divided into two categories, specimen geometry and 
material properties. In the specimen geometry category the two variables of interest are the taper ratio and 
aspect ratio. For convenience these definitions are repeated here, namely, 

W 

T.R. = — - (5- la) 

W, 

and 

A.R. = — - (5-lb) 

w 2 

Figure 5.1 is a schematic diagram which shows the test specimen geometry. As in Fig. 1.1, the dimensions 
W] and W 2 are the unsupported plate width at the narrow and wide (top and bottom) ends, respectively, and 
L is the unsupported plate length. The term unsupported is used because in order to provide the appropriate 
boundary conditions for the experiments, the specimens are oversize to accommodate fixturing. The gray 
areas, with associated dimensions, in Fig. 5.1 represent areas of the specimen held within the test fixture. 
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The dimensions Wj * W 2 , and L are the dimensions of the test specimens required to obtain the necessary 
fixturing areas and unsupported lengths. 



Figure 5.1 Test Specimen Geometry 

Table 5.1 presents the specific aspect ratio and taper ratio combinations investigated and the corresponding 
values for dimensions W j , W 2 , and L . 
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Table 5.1 Range of Geometric Parameters Studied and Specimen Dimensions 


■ Taper Ratio 

|iAspcctRatto i 

W, (mm) 

W 2 (nun) 


1.0 

1.0 

368 

368 

375 

1.5 

1.0 

247 

372 

375 

1.5 

2.0 

248 

370 

730 

3.0 

1.0 

125.6 

375 

375 


The effects that material fiber orientation angle and plate stacking sequence have on the nondimensional 
parameters were discussed in connection with Fig. 3.1. To investigate these effects, specimens for the 
experiments have ply orientations ranging from 0 to 90 degrees at 0, 30, 45, 60, and 90 degrees. To assess 
the effect of and D 26 , both clustered and interspersed stacking sequences are studied for each value of 8. 
As limiting cases to bound the problem, the [0] 2 4 and [90] 24 plates are studied. And finally, a quasi-isotropic 
laminate [445/0/-45/90] 3s is studied because it is a laminate of interest. Aluminum specimens are also tested. 
Table 5.2 presents the matrix of the specific ply orientation and stacking sequence combinations 
investigated. 


Table 5.2 Laminates Investigated 


0 (deg.) 

S kin Seq ' 

0 

[0] 2 4 

30 

[+30/-30] 6s 

30 

[+306/-30 6 ] s 

45 

[+45M5] 6s 

60 

[+60/-60] 6s 

60 

[+606/-60 6 ] s 

90 

[90] 24 

Quasi- 

Isotropic 

[+45/0/-45/90] 3s 

Isotropic 

aluminum 
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5.3 Specimen Fabrication 

Test panels are laid up by hand from Hercules AS4/3502 graphite/epoxy prepreg and cured in an autoclave 
using the manufacturers recommended cure cycle. A single test specimen is machined from each test panel 
using a diamond abrasive cutting wheel with flood coolant. Special care is taken to ensure that the ends 
remain as flat, parallel, and undamaged as possible. The Micro-Measurements CEA-125UR-350 [11] strain 
gages are applied using M-Bond 200 strain gage adhesive [12]. To enhance the shadow moire images, the 
front of the specimen is painted white. The sections outside of the unsupported area are masked off before 
painting to ensure that the test fixture contacts only unpainted surfaces. Prior to testing, thickness 
measurements are made at three locations along each of the four edges of the specimen. The average 
specimen thickness for each edge, and the resulting average ply thickness for the specimen are reported in 
Table 5.3. Recall from Table 2.2 that the analysis assumed a ply thickness of . 127 mm. 


Table 53 Average Specimen Thickness 



ip 


pftf 

1 Jt 

k-< 

iWiil 13 


BUSH 


1.0 

1.0 

[0] 2 4 

2.95 

3.00 

2.97 

2.97 

.124 

1.0 

1.0 

[90] 24 

2.79 

2.84 

2.87 

2.87 

.119 

1.0 

1.0 

[+30/-30] 6s 

2.87 

2.90 

2.90 

2.87 

.119 

1.0 

1.0 

[+30^-306], 

2.87 

2.84 

2.90 

2.87 

.119 

1.0 

1.0 

[+45/-45] 6s 

2.90 

2.92 

2.90 

2.90 

.122 

1.0 

1.0 

[+45/0/-45/90] 3s 

2.95 

2.90 

2.87 

2.92 

.122 

1.0 

1.0 

[+60/-60] 6s 

2.95 

2.95 

2.97 

2.97 

.122 

1.0 

1.0 

[+6<V-60 6 ] s 

2.84 

2.90 

2.84 

2.87 

.119 

1.0 

1.5 

[0] 2 4 

2.79 

2.87 

2.84 

2.84 

.119 

1.0 

1.5 

[901 24 

2.82 


2.84 

2.84 

.119 

1.0 

1.5 

[+30/-30] & 

2.90 

2.90 

2.90 

2.92 

.122 

1.0 

1.5 

[+ 30 ^- 30 ^ 

2.82 

2.87 I 

2.87 

2.84 

.119 

1.0 

1.5 

(+45/-45]6g 

2.92 

2.92 

2.90 

2.92 

.122 

1.0 

1.5 

[+45/0/-45/90] 3s 

3.00 

2.90 

2.87 

2.90 

.122 

1.0 

1.5 

[+60/-60] 6s 

.2.87 

2.84 

2.87 

2.90 

.119 
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Table 5.3 Average Specimen Thickness 


AJR. 

TJL 

Sequence 

ill f : f 

Top 

(mm*) 

Right 

(mm) 

r; 

• , ’. • « >!• 

thicknes 

s 

1.0 

3.0 

[0] 2 4 

2.87 

2.90 

2.90 

2.87 

.119 

1.0 

3.0 

[+30/-30] 6s 

2.87 

2.90 

2.90 

2.90 

.119 

1.0 

3.0 

[+306/-306]* 

2.84 

2.84 

2.79 


.117 

1.0 

3.0 

[+45/-45] 6s 

2.87 

2.90 

2.90 

2.90 

.122 

1.0 

3.0 

[+45/0/-45/90] 3s 

3.02 

3.05 

3.05 

3.02 

.127 

1.0 

3.0 

[+6(y-60 6 ] s 

2.90 

2.92 

2.92 

2.95 

.122 

2.0 

1.5 

10124 

2.87 

2.90 

2.87 

2.92 

.122 

2.0 

1.5 

[90324 

2.92 

2.95 

2.90 

2*95 

.122 

2.0 

1.5 

[+30/-30] 6s 

2.95 

2.95 

2.95 

2.97 

.122 

2.0 

1.5 Vj 

[+306/-30 6 ] 8 

2.64 

2.67 

2.62 

2.62 

.109 

2.0 

1.5 

[+45/-45J6S 

x- 2.92 

2.92 

2.87 

2.90 

.122 

2.0 

1.5 

[+45/QM5/90] 3s 

3.18 

3.23 

3.20 

3.23 

.135 

2.0 

1.5 

[+60/-601& 

2.97 

2.95 

2.92 

2.95 

.122 

2.0 

L5 

[4^60^ 

2.79 

2.84 

2.74 

2.79 

.117 


5.4 Test Fixture 

Specialized test fixtures are used to provide the required boundary conditions on the wide variety of tapered 
specimens, and to ensure uniform loading of the specimen. The main components of the test fixture are: top 
and bottom clamping blocks with sliding plates to enforce the clamped boundary at the top and bottom of 
the plate, knife edge rails to provide simple supports along the sides of the plate, knife edge restraints to 
keep the knife edge rails from moving out-of-plane during the test, and a top loading platen loaded through a 
semi-spherical ball to provide uniform loading by accommodating any lack of parallelism of the ends of the 
specimens. Figure 5.2 is an isometric drawing showing the major components of the test fixture and a cross 
section of the knife edge rails. 
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Figure 5.2 Schematic of Test Fixture 

Figure 5.3 is a photograph from the front of the test fixture with a specimen. The taper ratio of the specimen 
is 1 .5 and the aspect ratio is 1 .0. Visible are the top and bottom clamping blocks, the right and left knife edge 
rails, top and bottom loading platens, and the right and left DCDT mounts, to be discussed shortly. The top 
loading platen and the semi-spherical loading ball, which will also be discussed shortly, are visible as well. 
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The specimen is painted white to enhance the shadow moire patterns, and strain gage lead wires running 
vertically are visible on the face of the specimen. 



Hemispherical Ball 
Top Loading Platen 


DCDT Mount 

Top Clamping 
Block 

Knife Edge Rail 


Bottom Clamping 
Block 

Bottom Loading 
Platen 


Figure 5.3 Front View of Test Fixture 
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Figure 5.4 is a close-up photograph showing the lower left-hand comer of the test fixture. Visible are the 
interactions of the bottom clamping block and sliding plate, knife edge rail, and knife edge restraints. 



Knife Edge Rail 


— Knife Edge Restraint 


— Bottom Clamping 
Block 

— Bottom Sliding Plate 


— Knife Edge Restraint 


Figure 5.4 Detail of Test Fixture 


The bottom clamping mechanism consists of the bottom clamping block and a bottom sliding plate. The 
bottom 9.53 mm of the specimen is clamped between the clamping block and the sliding plate. The sliding 
plate is in the foreground of Fig. 5.4 and the gap due to the specimen thickness is visible as a dark vertical 
line. The gap between the knife edge rail and bottom clamping block/bottom sliding plate assembly due to 
the 3.18 mm spacer is also visible. Before testing the spacer is removed and this gap allows for end- 
shortening of the specimen during application of the compressive load. The knife edge rails grip the outer 
6.35 mm of the specimen, and one knife edge restraint bolts to the clamping block while the other bolts to 
the sliding plate. As can be seen in Fig. 5.2, the specimen is squeezed between knife edges on the knife edge 
rail and the sliding bar. The interaction between fixture components shown in Fig. 5.4 occurs at all four 
comers of the specimen. 
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Figure 5.5 shows the detail of the semi-spherical ball arrangement used to allow the loading platen to 
accommodate any lack of parallelism between the top and bottom edges of the specimen 



Top Crosshead of 
Load Frame 

Rectangular 

Block 

Dowel Pin 
(1 of 4) 

Convex Half of 
Semi-Spherical 
Ball 


Concave Half of 

Semi-Spherical 

Ball 

Top Loading 
Platen 


Figure 5.5 Detail of Semi-Spherical Loading Ball 

The rectangular block is indexed in the center of the load frame crosshead with small blocks that lit into a 
square opening normally used for tensile wedge gups. The convex half of the semi-spherical ball is aligned 
in the center of this block with four dowel pins. The concave half of the ball fixture is indexed to the center 
of the 38.1 mm thick top loading platen by two other dowel pins (not visible in Fig. 5.5). The four long bo 
visible in the upper portion of Fig. 5.3 are designed to hold the loading platen up when the specimen is 
unloaded. The entire test fixture is centered front ,0 back in the load frame wifi, mechanical stops on the top 
and bottom loading platens. All of this alignment is ,0 ensure that the load is transmitted from the center of 
the crosshead through the center of the top loading platen, through the center of the test fixture, and 
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therefore, through the center of the specimen. This alignment is critical. The front to back alignment is 
especially critical for preventing bending, an event which would seriously affect the buckling loads. 

5.5 Instrumentation 

A variety of sensors are used to measure the response of the plates. Load is measured using the standard load 
cells of the Tmius Olsen load frame. Specimen end-shortening is measured, relative to the crosshead support 
columns of the load frame, on both sides of the specimen with direct current displacement transducers 
(DCDTs). End-shortening measurements on both the left and right sides of the specimen reveal any side to 
side nonuniformities in the loading of the specimen. The barrels of the DCDT’s clamp to the support 
columns with blocks that are adjusted up/down and in/out to accommodate specimens of varying length. The 
cores of the DCDTs attach to the top loading platen with swivel joints that allow the core to remain vertical 
even if the loading platen should rotate slightly due to a specimen with non-parallel ends. Coarse zero 
adjustment for each DCDT is provided by sliding the clamping block up or down on the support column, 
and fine zero adjustment is provided by turning the threaded rod connecting the core to the loading platen. 
Figure 5.6 shows the detail of the DCDT mount. 
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Swivel Joint 


Top Loading Plate 


Figure 5.6 Detail of DCDT Mount 

Out-of-plane deflection of the specimen is measured by a third DCDT. This DCDT is oriented horizontally 
and is spring loaded to maintain contact with the center of the plate. Inplane strains are measured with 350 
ohm foil-gage strain gages adhesively bonded to the specimen. Back-to-back strain measurements are taken 
at all four corners and the center of the plate. The strain measurements at the comers of the plate are useful 
to determine if any undesirable bending due to nonuniform loading occurs prior to buckling. The strain 
measurements at the center of the plate provide additional data to help determine the buckling load. All data 
is acquired using an Orion Solartron data logger controlled by an IBM PC. Shadow moir6 is used to 
visualize the buckled mode shapes. The specimen grating has a frequency of 1.968 lines/mm and the lines 
are aligned vertically. The angle of incidence of the light is 45 degrees, and the angle of observation is zero 
degrees. A Sony video camera and U-Matic video recorder are used to record this video data. Figure 5.7 
shows the location of the 10 strain gages and the out-of-plane DCDT. The rectangles signify the location of 
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the strain gages. The odd numbered strain gages are on the front (painted white) of the specimen, while the 
even numbered strain gages are in the same location on the back (unpainted) of the specimen. The solid 
circle represents the location of the out-of-plane DCDT. 



Figure 5.7 Shadow Moire Setup, Strain Gage and Out-of-Plane DCDT Locations 

5,6 Experimental Procedure 

The following is the list of steps used to load the specimen in the test fixture, ensure alignment and 
uniformity of loading, and test the specimen: 

1. Clamp the bottom 9.53 mm of the specimen loosely in the bottom clamping block and sliding plate 
mechanism using a bar clamp near the left and right edges of the specimen. (A bar clamp is similar to 
a ‘C’ clamp.) 
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2. Adjust the specimen until the distance from the bottom edge of the specimen to the edge of the bottom 
clamping block is the same, to within 0.25 mm, on both the right and left edges of the specimen. This 
measurement is made by using a steel rule. 

3. Tighten the bar clamps to firmly grip the specimen. Tighten the bolts holding the bottom sliding plate 
to the bottom clamping block. 

4. Place the 3.18 mm spacing pads next to both the left and right bottom edges of the specimen. 

5. Rest the bottoms of the knife edge rails on the spacing pads and clamp them loosely to the edges of 
the specimen using the sliding bar shown in Fig. 5.2. Inspect to make sure that the specimen is fully 
inserted into the knife edge rails. 

6. Firmly clamp the bottom knife edge restraints against the knife edge rails using the bar clamps. 

7. Tighten all of the bolts on the knife edge rails and the bolts that secure the bottom knife edge restraints 
to the bottom clamping block and bottom sliding plate. 

8. Place the top clamping block and sliding plate assembly on the top of the specimen and loosely clamp 
the specimen in the assembly using the bar clamps. Measure the space between the top specimen edge 
and the top clamping block edge on both the left and right edges of the specimen, and adjust until 
there is less than 0.25 mm difference between the right and left edges, as was done for the bottom 
clamping block. Firmly clamp the specimen, and tighten the bolts holding the top sliding plate to the 
top clamping block. 

9. Clamp the top knife edge restraints against the knife edge rails and tighten the bolts that secure the top 
knife edge restraints to the top clamping block and top sliding plate. 

10. Position the fully assembled fixture against the alignment stops on the top and bottom loading platens. 

11. Align the marked center of the fixture with the marked centerline of the bottom loading platen. 

12. Load the specimen to 1334 Newtons. 

13. Loosen and then retighten finger tight the four nuts on the bolts holding up the top loading platen. 

14. Ensure that the knife edge rails are perpendicular (front to back) with the top and bottom clamping 
blocks. 

15. Unload the specimen. 
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16. Attach all of the strain gage leadwires to the data acquisition system. 

17. Start the data acquisition system, slowly load the specimen to 2220 Newtons, then unload the 
specimen. 

18. Examine the response of the left and right DCDT’s for uniformity of end-shortening from left to right 
on the specimen (This checks for misalignment of the top loading platen). 

19. Examine the response of the four strain gages mounted at the top of the plate. Again check for 
uniformity of loading right to left. Also check for front to back bending. 

20. Make adjustments to fixture and specimen, if necessary, to attain uniformity of loading. 

21. Place the shadow moir6 grating in front of the specimen 

22. Start the video recorder and data acquisition system. 

23. Slowly load the specimen to a load level approximately 30% greater than the load at which a 
significant number of shadow moi r€ fringes first start to appear. 

24. Unload the specimen. 

5.7 Determining Buckling Load 

Buckling is characterized experimentally by a change in slope of the load vs. end-shortening curve. This 
change of slope of the load vs. end-shortening curve occurs after the plate buckles because the out-of-plane 
deflection of the plate has a softening effect on the structure. This results in significantly more end- 
shortening for a given increase in load. The experimental buckling load is determined by the intersection 
point of straight lines fit through the prebuckling and postbuckling portions of the load vs. end shortening 
relations. Figure 5.8 presents typical load vs. end-shortening relations as measured by the left and right 
DCDTs. The particular specimen for which the data are presented has aspect ratio 1 . 0 , taper ratio 1 . 5 , and 
stacking sequence [+60£/-60 6 ] s . The straight lines through the prebuckling and postbuckling portions of the 
relations shown in Fig. 5.8 are faired in by hand. While this technique seems quite arbitrary, using the 
computer to calculate a linear least-squares fit would also require arbitrary choices for locations to begin and 
end the fits to both the prebuckling and postbuckling portions of the relations. Drawing the lines by hand 
also allows for judgement to avoid the initial portions of the curves when the fixtures are taking up slack. 
Figures 5.9 - 5.12 present load vs. strain and load vs. out-of-plane relations for the same specimen. 
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Strain Ggage #2 
Strain Ggage #3 
Strain Ggage #4 


Strain (|i strain) 

Figure 5.9 Typical Load vs. Strain Relations: Top Corner Gages 
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Strain (| l strain) 

Figure 5.11 Typical Load vs. Strain Relations: Center Gages 
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Figure 5.12 Typical Load vs. Out-of-Plane Displacement Relation, Center of Specimen 

Note that in Fig. 5.8 there is very little difference between the load vs. end shortening relations obtained 
from the left and right DCDTs, indicating that the specimen was loaded uniformly across its width. The load 
vs. strain relations in Figs. 5.9 -5.10 also indicate a uniform load across the specimen width, in addition to 
revealing minimal bending of the specimen before buckling. Finally, note that the load at which the back-to- 
back strain gages in Figs. 5.9 - 5.11 begin to diverge and the load at which the out-of-plane displacement 
begins to have a large increase are approximately the same as the buckling load determined using the 
outlined procedure in connection with Fig. 5.8. 

It must be stated that the process of fitting a line through the prebuckling and postbuckling portions of the 
load vs. end-shortening relationships is not always as clear as the case shown in Fig. 5.8. For the stiffer 
laminates the change in slope between the prebuckling and postbuckling stiffness is not as great. Figures 
5.13 and 5.14 demonstrate this point by presenting load vs. end shortening curves for a [±30] 6s laminate and 
a [0]24 laminate, respectively. Comparing Figs. 5.8 and 5.13 reveals that the difference between prebuckling 
and postbuckling slopes is less for the [±30]5 S laminate. Comparing Fig. 5.8 with Fig. 5.14 shows that this 
trend continues for the [0] 24 laminate. Note also that the uncertainty is amplified further for this laminate 
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because the total end shortening is much smaller due to its high stiffness parallel to the loading direction 
compared to both the 30 degree and 60 degree laminates. 



Figure 5.13 Load vs. End-Shortening Relations: [±30] 6s Laminate 
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5.8 Calculating Experimental Buckling Coefficient 

Due to the definition of the buckling coefficient in Ch.3, there is a subtle difference in calculating the 
buckling coefficient for experimental results compared to calculating the buckling coefficient for the 
analytical results. As shown in Eqs. 2.128 *2.131, the eigenvalue computed by the semi-analytic approach is 
the loading parameter p which has the form 


where 


tt 7 d„D 22 


N‘ 



(5.2) 


(5.3) 
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Remember that Wj and bj are the unsupported plate width and plate half- width, respectively, at the narrow 
end of the plate. In order to pose the analytical results in the form of Eq. 3.5, p is modified according to 


N 


P = 
4p = 

4p = 
W 


2 - 


w 


This yields 


n 7^11^22 


K w ? 


analym 1 


7i 2 7D u D 22 

(5.4) 

PWj 


n jD u D n 


PW 2 


7i ^/D n D 22 


. - if«=0 ■ 

(5.5) 


When computing the experimental buckling coefficient, however, the plate width at the narrow end of the 
plate has dimension Wj , which is wider than the unsupported plate width Wj, as shown in Fig. 5.1. This 
changes the above derivation according to 


n V D 11 D 22 

n; w? 

_ txp 

^ 2 V D U D 22 

PW* 

W,* 2 Ad 22 



PW 

1^ - 2 1 

(?) 

lw, J 

I* 2 

71 j 

'd„d 22 


(5.6) 
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to 


w 2 / w r 

This shows that the measured buckling loads must be multiplied by the quantity . =- 

provide correct buckling coefficients for comparison with the buckling coefficients calculated by the semi- 


analytic approach. 


The next chapter describes experimental results and offers a comparison with semi-analytical predictions. 
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Chapter 6 

Experimental Results and Comparison with Analysis 


In this chapter, comparisons between the measured and predicted prebuckling stiffnesses are made for the 
ranges of plates tested. Buckling loads are compared, as are buckled mode shapes. As will be seen, the 
agreement between predictions and observation is generally quite good. However, for selected cases the 
agreement is not quite as good. The final section of this chapter offers possible explanations for the 
discrepancies that do occur. 

6.1 Comparison of Experimental and Analytical Prebuckling Plate Stiffness 

The formula for prebuckling plate stiffness, i.e., the slope of the prebuckling load vs. end-shortening 
relationship, was developed in Ch. 2 based on the assumed prebuckling inplane force resultant distribution, 
and is given by Eq. 2.141. Table 6.1 compares the calculated and experimentally determined prebuckling 
plate stiffness values for all of the plate geometry and stacking sequence combinations tested. The 
experimental stiffness values reported are the slopes of the lines fit by hand through the prebuckling portion 
of the load vs. end-shortening relationships, as described in Ch. 5. With the exception of the [0] 24 laminates, 
there is very good agreement between the calculated and experimental values. To a large degree the table 
speaks for itself. Typically the calculated values and experimental results are within 9%. However, there is 
between 27% and 35% difference for the [0] 24 laminates. Furthermore, the differences for all cases 
considered seem to be random in that the calculated results are neither always greater nor always less than 
the experimental results. 
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Table 6.1 Calculated and Experimental Prebuckling Plate Stiffness Values 


Sequence v : 

lii&i 

If* 

m 

lllr 

Experimental 

mz:.smmms 

. :->« . .X X - x . . X ; X- , ■ j. .• 

(N/ram) 

1 • | 

[0] 2 4 

1.0 

1.0 

3.89E5 

3.05E5 

[0124 

1.5 

1.0 

3.21E5 

2.5 1E5 

[0] 24 

3.0 

1.0 

2.37E5 

1.75E5 

{+30/-30]6s 

ill 

1.0 

1.50E5 

1.56E5 

[+30/-301& 

iiii 

l!p|#|P| 

frf#?1.24E5 : : 

1.23E5 

[+30/-301& 

sliill 

iilll 

?: I3E4 

9.02E4 

[+30/-30J& 

llill! 

•lA^Olp 

|jgjili8E4 

6.58E4 

[+30^-30*], 

1.0 

1.0 

1.50E5 

1.56E5 

[ + 30g/-30 6 ] s 

1.5 

1.0 

1.24E5 

1.24E5 

t + 30 6 /-30 6 ] s 

3.0 

1.0 

9.13E4 

9.34E4 

[ + 30 6 /-30 6 ] s 

1.5 

2.0 

6.18E4 

6.31E4 

Jf445M5] & 

10 

1.0 

‘i. " '■ 

6.06E4 

6.59E4 

r+4S/-45] 6s 

f iii 

1.0 

, wwmm:.:*. ■ 

4.98E4 

mrnmmmta ... : 

5.57E4 

t+45/-45] 6s 

ilii 


3.68E4 

3.96E4 

| 1+45A45J& 

1.5 

2.0 

2.49E4 

2.60E4 
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9.39E4 
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1.69E4 

[+60 6 /-60 6 ] s 

1.0 

1.0 

3.84E4 

3.85E4 

[+60 6 /-60 6 ] s 

1.5 

1.0 

3.15E4 

3.15E4 

[+606/-60 6 ] s 

3.0 

1.0 

2.33E4 

2.29E4 

[+606/-606], 

1.5 

2.0 

1.57E4 

1.79E4 
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Table 6.1 Calculated and Experimental Prebuckling Plate Stiffness Values 


IS 



1.5 1.0 2.73E4 2.56E4 

190124 

1.5 2.0 1.38E4 1.37E4 


6.2 Comparison of Experimental and Analytical Buckling Loads 

As mentioned in Ch. 5, several aluminum plates were tested to provide a baseline comparison. The buckling 
coefficients from the experiments with aluminum plates are given in Table 6.2 along with the buckling 
coefficients for the aluminum plates as computed by the semi-analytical analysis. In Table 6.2, and in similar 
tables to follow, finite element results are included because not all of the experimental cases are considered 
in the figures of Ch. 4. In this way, future comparisons between the semi-analytic approach and the finite 
element results are available. Discussion in this chapter, however, will be limited to comparing experimental 
results to results from the semi-analytic analysis. Additionally in Table 6.2, the number of halfwaves in the 
buckled mode shape is given by numbers in parentheses in the Experimental and Analysis columns. The 
fundamental case for a baseline comparison is a square (A.R. = 1.0, T.R. =0) isotropic (aluminum) plate 
because it simplifies all of the material and geometric parameters in the analysis. Table 6.2 shows 
approximately 1.5% disagreement for this important case. The other two rectangular cases (A.R. = 1 .0) have 
disagreements of 10% and 2.2%. As the taper ratio increases, the effect of the geometric parameters in the 
analysis is investigated, while the material parameters are still simplified. The moderately tapered geometry 
(A.R. = 1.0, T.R. = 1.5) shows a 4.1% disagreement and the extremely tapered geometry (A.R. = 1.0, 

T.R. = 1.5) shows a 1.2% disagreement. Nothing can be stated regarding the sign of these disagreements, in 
some cases the analytical results are conservative and in other cases they are nonconservative. In all cases 
the number of half waves observed is correctly predicted. 
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Table 6.2 Experimental Buckling Loads for Aluminum Plates 



gg 

u 

gil 

Pi! ■ i 

r Jpif 

Aluminum 

1.0 

1.0 

6.6(1)* 

6.7(1) 

6.8 

Aluminum 

1.0 

1.5 

4.9(1) 

5.4(1) 

5.4 

Aluminum 

1.0 

3.0 

4.5(3) 

4.4(3) 

4.5 

Aluminum 

1.5 

1.0 

7.3(1) 

7.0(1) 

7.0 

Aluminum 

3.0 

1.0 

7.8(1) 

7.7(1) 

8.2 


* The number of halfwaves in the budded mode shape is indicated by the number in parentheses. 


The next logical step, after investigating the baseline aluminum material, is to investigate the effects of the 
orthotropic parameters independently from the effects of the anisotropic parameters. Table 6.3 presents 
analytical and experimental results for [0] 24 and [90] 2 4 composite plates. As mentioned in Ch. 3, these 
laminates are limiting cases for the orthotropic parameters in this analysis. As can be seen from Table 6.3, 
there is much better agreement, in general, between the calculated and measured buckling coefficients for 
the [90] 24 composite plates than for the [0] 24 plates. Some difficulties due to initial imperfections were 
encountered while conducting the experiments on the [0] 2 4 laminates. The disagreements between the 
experimental and analytical results range from 22% to 32% for the [0] 24 laminate, compared to a range from 
0% to 7.8% for the [90] 24 laminate. 

Unidirectional composites have inherent difficulties as test specimens due to the large difference in stiffness 
between the two principal directions. Unidirectional plates can be less flat than angle ply or cross ply 
laminates due to a high sensitivity to fiber misalignment as the individual plies are laid down during 
fabrication. Indeed, one of the [0] 24 specimens was completely discarded because the specimen showed 
significant shadow moire fringe patterns in the unloaded state and exhibited large out-of-plane 
displacements immediately upon application of the load. Also, not having the 0 degree direction perfectly 
perpendicular to the loaded edges can be a problem, i.e., the laminate could actually be a [3] 24 instead of a 
[0] 24 laminate. Additionally, it would appear that the difficulties associated with testing unidirectional 
specimens are more severe when the specimens are loaded in the stiff direction than when they are loaded 
perpendicular to the stiff direction. Recall from Ch. 2 that the assumed prebuckling inplane stress 
distribution is more accurate for the specimens loaded perpendicular to the stiff direction. Note the 
disagreement in the number of halfwaves predicted by the analysis for the long [90] 24 laminate. 
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Table 6.3 Experimental Buckling Loads for Unidirectional Composite Laminates 


i m 
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|||§ 

wm 

III II 


m 

[0] 2 4 

1.0 

1.0 

11.0(1) 

14.6(1) 

12.7 

[0124 

3.0 

1.0 

8.7(1) 


9.2 

H^jlP 

mat 

1.0 

3.8(2) 

mm& 

3.9 

■B 

mil 

1.0 


4.2(2) 


{SO}* 

KB 

2.0 

3.8 (5) 

■BE m 

3.5 


The final step in the investigation is including all of the effects of the geometric, orthotropic, and anisotropic 
parameters. To that end, buckling coefficients for angle-ply composite plates are presented in Table 6.4. As 
in Table 6.3, there is better agreement between the measured and calculated buckling coefficients for 
specimens loaded perpendicular to their stiff direction. The [+30/-30]g and [+30g/-306] s have differences up 
to 32% between the measured and calculated buckling coefficients, while the largest difference for the [+60/ 
-60] 6 and [+606/-606],; laminates is 13%. For 75% of the cases shown in Table 6.3 there is less than 12% 
difference between the measured and calculated buckling coefficients. In only 4 of the 19 cases shown is the 
difference greater than 20%. Note that for the long plates, A.R. = 2.0, the analysis underpredicts the number 
of halfwaves by 1 . It is also worth noting, however, that the cases in which there is greatest disagreement 
between analytical and experimental buckling coefficients are not always those in which there is 
disagreement in the number of halfwaves in the buckled mode shape. 


Table 6.4 Experimental Buckling Loads for Angle-Ply Laminates 
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Table 6.4 Experimental Buckling Loads for Angle-Ply Laminates 


Stacking 

- v ' 

: m 

Experimental 

■ m - ! 

‘kniyiisf: 

wmm 

. xi:- xix :■ 'x x ■ • • x* •• 

[±45] 6s 

1.0 

1.0 

8.7(1) 

9.2(1) 

9.2 

[±45] 6s 

1.5 

1.0 

9.0(1) 

9-8(1) 

10.0 

[±45] 6s 

3.0 

1.0 

11.4 (2) 

10.6(1) 

12.1 

[±45] 6s 

1.5 

2.0 

8.9 (3) 

9.3 (2) 

8.9 

■■■ 

1.0 

IIP 

7.9(2) 

8.4(1) 

8.3 

{±60]^ 

| llil 

IIP 

8.6(1) 

8.9(1) 

illPli 

[±60165 

1.5 

2.0 

7.8(4) 

8.8(3) 

8.7 

[+606/-60 6 ] s 

1.0 

1.0 

6.4(1) 

6.9(1) 

6.8 

[+606/-60 6 ] s 

1.5 

1.0 

6.0 (2) 

6.8(1) 

6.8 

[+606/-60 6 1 s 

3.0 

1.0 

6.6 (2) 

7.5 (2) 

7.7 

[+606/-60 6 ] s 

1.5 

2.0 

5.7 (4) 

5.8(3) 

5.7 


Table 6.5 presents analytical and experimental results for quasi-isotropic composite plates. The 
disagreements associated with these laminates are comparable to those of the [±45]6 S laminates shown in 
Table 6.4. In all cases the number of half waves predicted was observed. 


Table 6.5 Experimental Buckling Loads for Quasi-Isotropic Laminates 
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Stacking 

T.R. 

p| | 

■ ' I . 

Experimental 

' kc : 

gg| 

Si: 

[+45/0/-45/90] 3s 

1.0 

1.0 

7.1(1) 

7.4(1) 

7.4 

+45/0/-45/90] 3s 

1.5 

1.0 

8.2(1) 

7.7(1) 

7.5 

+45/0/-45/90] 3s 

3.0 

1.0 

9.6(1) 

8.1(1) 

9.9 

+45/0/-45/90] 3s 

1.5 

2.0 

7.1 (2) 

6.3 (2) 

6.3 


6.3 Comparison of Experimental and Analytical Buckled Mode Shapes 

This section presents representative fringe patterns to demonstrate the effects of plate geometry and ply 
stacking sequence on the buckled mode shape. The drastic effect the direction of orientation of the principal 
direction having the greater stiffness has for laminates that are highly orthotropic is shown in Fig. 6.1. The 
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specimen on the left has a [ 0]24 stacking sequence while the specimen on the right has a [90]24 stacking 
sequence. Both specimens have aspect ratio 2.0 and aspect ratio 1.5. Notice that because the [90] 2 4 laminate 
is relatively soft in the direction of the applied load, it buckled into 5 half-waves along the length. The [0]24 
specimen, on the other hand, is very stiff in the direction of the applied load so it buckled into only a single 
half wave. Because neither of these specimens has a D ^ or D 26 component in the bending stiffness matrix, 
the buckled mode shapes are not skewed. 



Figure 6.1 Experimental Buckled Mode Shapes for [0]24 (left) and [90)24 (right) Laminates: 

A.R. = 2.0, T.R. = 1.5 
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The effect of the D !6 and D 26 terms skewing the buckled mode shape is visible in Fig. 6.2. The specimen on 
the left has stacking sequence [+60/-60] 6s and the specimen on the right has stacking sequence [+60 6 y-60 6 ] s . 
The specimens both have aspect ratio 2.0 and aspect ratio 1.5. Note that both specimens buckled into the 
same number of half-waves, but the waves are not symmetric from left to right for the specimen on the right. 
As was seen in Fig. 3.1, the [+605/-606] s laminate has much larger values for D 16 and D 26 . However, D 16 
and D 2 6 are not zero for the [±60]g s laminate and close inspection of the left portion of Fig. 6.2 shows very 
mild skewing for this specimen also. 



Figure 6.2 Experimental Buckled Mode Shapes for [±60]^ (left) and [+60 6 /-60 6 ] s (right) Laminates: 

A.R. = 2.0, T.R. = 1.5 
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Figure 6.3 shows a similar effect for [+30/-30] 6s and [+30g/-306] s laminates having aspect ratio 2.0 and 
aspect ratio 1.5. Note that for this case, the f+SO^-SOgJs laminate buckles into 3 half-waves while the [+30/- 
30 ] 6s laminate buckles into only 2 half-waves. The skewing for the t+SO^-SOgls laminate is quite extreme. 



Figure 6.3 Experimental Buckled Mode Shapes for [±30]^ (left) and [+30j/-30j] s (right) Laminates: 

A.R. = 2.0, T.R. = 1.5 
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Figure 6.4 shows a specimen with stacking sequence [ +30(/-30(J s , but with plate taper ratio 3.0 and aspect 
ratio 1.0. It is important to note the large effect the combination of an extremely tapered geometry and a 
highly anisotropic material have on the buckled mode shape. There is clearly a high lack of symmetry for 
this case. 



Figure 6.4 Experimental Buckled Mode Shape [+30$/-30£] s Laminate: A.R. = 1.0, T.R. — 3.0 

Although the number of halfwaves calculated by the analysis are presented in Tables 6. 1 through 6.4, for 
comparison purposes Figs. 6.5 through 6.8 present contour plots of the buckled mode shape calculated by 
the semi-analytic analysis for the same plate geometries and stacking sequences shown in Figs. 6. 1 through 
6.4. Comparing the two sets of mode shapes reinforces what was shown in Tables 6.1 through 6.4 regarding 
the number halfwaves in the buckled mode shape computed by the analysis compared to the number of 
halfwaves exhibited by the specimens during testing. Specifically, for the longer aspect ratios the analysis 
underpredicts the number of halfwaves in the fundamental buckling mode. What is not apparent from the 
tables, however, is that the amount of skewing in the calculated buckled mode shapes is very similar to that 
shown in the experiments. 


6.3 Comparison of Experimental and Analytical Buckled Mode Shapes 


127 



Chapter 6: Experimental Results and Comparison with Analysis 



Figure 6.5 Calculated Buckled Mode Shapes of [O]^ (left) and [90]24 (right) Laminates: 
A.R. = 2.0, T.R. = 1.5 (compare with Fig. 6.1) 
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Figure 6.6 Calculated Buckled Mode Shapes for [±60] 6s (left) and [+60 6 /-60 6 ] s (right) Laminates: 

A.R. = 2.0, T.R. = 1.5 (compare with Fig. 6.2) 
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Figure 6.8 Calculated Buckled Mode Shape [+30</-306] s Laminate: A.R. — 1.0, T.R. — 3.0 

(compare with Fig. 6.4) 


6.4 Summary of Sources of Disagreement between Experiments and Analysis 

Many possible sources of disagreement have been discussed throughout this chapter and throughout this 
investigation. The ones believed to have the greatest influence on the experimental and analytical results are 
summarized in this section. 

Recall from Ch. 2 the discussion associated with Figs. 2.6 - 2.17 regarding the assumed form of the 
prebuckling inplane force resultant distribution. In that discussion it was pointed out that the discrepancy 
between the assumed form and the finite element results appeared to be a strong function of the degree of 
orthotropy and the degree of taper. This discrepancy in prebuckling force resultant distribution appears to be 
a major contributor to the disagreement. Recall the isotropic cases show much better agreement than the 
composite laminates, and the composite laminates that have their stiff axis oriented perpendicular to the 
loading direction show better agreement than the composite laminates that have their stiff axis oriented 
parallel to the direction of loading. This dependence of the orientation of the stiff axis was also evident when 
comparing buckling predictions from the semi-analytic analysis with buckling predictions from the finite 
element analysis. It can be expected, then, that for those cases where the prebuckling assumptions in error, 
differences between experiments and analysis will occur. 


6.4 Summary of Sources of Disagreement between Experiments and Analysis 


131 



Chapter 6: Experimental Results and Comparison with Analysis 


Recall from Ch. 5 that the straight lines through the prebuckling and postbuckling portions of the load vs. 
end-shortening relationships were fit by hand. For some experimental cases this was very straightforward 
while others required quite a bit more judgement. For those requiring more judgement, deviations between 
experiments and analysis can be expected. 

Recall from Ch. 2 that the material properties reported in Table 2.2 are typical values. Implicit in these 
values is a specified fiber volume fraction and resulting ply thickness. If the actual ply thickness is less than 
the specified ply thickness, the resulting laminate will have both a higher fiber volume fraction and a smaller 
plate thickness. The higher fiber volume fraction enhances the laminate stiffness and will contribute to a 
higher buckling load. The smaller plate thickness degrades the plate stiffness and will contribute to a lower 
buckling load. The opposite is true if the actual ply thickness is greater than the specified ply thickness. 

Any manufacturing errors involved in machining the ends of the specimens flat and parallel will manifest 
themselves in the experimental results. This is a much more important factor for the laminates that have their 
stiff direction oriented parallel to the direction of load. 

Finally, although the experimental fixtures appeared to do a good job providing uniform loading, there is no 
way the ensure a perfectly clamped boundary conditions on the ends or simply supported boundary 
condition on the edges. Unfortunately, is also no way to measure how well the fixture is providing the 
clamped or simply supported boundary conditions. 

In spite of difficulties and possible error sources, there is excellent agreement between the experiments and 
analysis regarding the prebuckling plate stiffness. There is also an acceptable level of agreement between the 
experiments and analysis to verify that the semi-analytic analysis is accurately predicting the buckling loads. 
The [0] 2 4 laminate, however, is the obvious exception to both of these statements. 
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7.1 Summary 

This research has used three different tools to investigate the budding response of symmetrically laminated 
composite plates having a planform area in the shape of an isosceles trapezoid. The nonparallel edges of the 
plate were assumed to have simply supported boundary conditions, while the parallel edges of the plate were 
assumed to be either simply supported or clamped pairs. The loading was assumed to be an inplane 
unidirectional force perpendicular to the parallel edges of the plate. First, a semi-analytic closed-form 
solution based on energy principles and the Trefftz stability criterion was derived in terms of six 
nondimensional parameters. Two of these parameters were measures of material orthotropy, two were 
measures of material anisotropy, and two were measures of plate geometry. Standard classical lamination 
theory, based on Kirchhoff s plate theory, was used in conjunction with the von Karman nonlinear strain- 
displacement relations. A simplified prebuckling force resultant distribution was assumed to simplify the 
analysis, and the out-of-plane displacement was approximated by a double trigonometric series. The 
Rayleigh-Ritz method of solution provided a symmetric eigenvalue problem that can be implemented on any 
computer. The second tool used in this investigation was a finite element analysis to compare buckling loads 
for several specific material and plate geometry combinations. Models consisting of a 10 by 10 mesh of 8- 
node shell elements were constructed using PATRAN and buckling loads were computed using the general 
purpose finite element code ABAQUS. The third tool used in this research was an experimental 
investigation. Special test fixtures were designed, and an experimental technique developed to verify the 
results obtained from the semi-analytic analysis. 

7.2 Conclusions Regarding the Analysis 

Intrinsic elements of the semi-analytic analysis were the inclusion of the effects of the anisotropic material 
terms (D 16 and D 2 6), the simplified prebuckling inplane force resultant distribution, and the series 
approximation for the out-of-plane displacement of the buckled mode shape. Including the anisotropic terms 
added complexity to the analysis because these terms were coefficients of mixed partial derivatives of the 
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out-of-plane displacement. Including these terms, however, was important because they were shown to 
decrease the buckling load, and cause the buckled mode shapes to be skewed. 

Rather than solving the first variation of total potential energy, the prebuckling force resultant distribution 

P 

was assumed to have the simplified form N = and N = N „ = 0, where P is the applied load 

x W(x) y xy 

and W(x) is the plate width. Comparison with finite element results showed that this was a very accurate 
approximation for isotropic material, and was more accurate for composite plates if the stiff material axis is 
oriented perpendicular to the loading direction compared to if the stiff material axis was oriented parallel to 
the loading direction. In general, the assumed force resultant values were within 1% of the finite element 
results, and in the worst case, the results were within 25%. No general trend existed as to whether the 
assumed force resultant distribution was conservative, or nonconservative. 

The series approximation for the out-of-plane displacement is more complex for tapered plates when 
compared to approximations for rectangular plates because the width of the plate is a function of position 
along the length of the plate. Including the anisotropic material terms added to this complexity, also, because 
the resulting skewed mode shapes prevented assuming any symmetry. The out-of-plane displacement was 
approximated by a double trigonometric series comprised of products of odd and even functions (i.e., 
symmetric and antisymmetric functions) that were assumed to have different forms for the two different sets 
of boundary conditions on the parallel edges of the plate. For the plate geometries considered, this series 
converged using 4 terms for isotropic materials. Including 5 terms in the series resulted in a 50 by 50 
eigenvalue problem that ran in less than 1/2 second on a Cray YMR Because of the fast convergence rate, 
and the fact that this analysis can run on any computer capable of supporting FORTRAN, this analysis can 
be a valuable design tool for the preliminary design of tapered composite plates. 

Based on comparison with finite element results, it appears safe to say that, for laminates of practical 
interest, this analysis can predict buckling coefficients accurate to within 11% for highly tapered plates and 
within 3% for moderately tapered plates. For more traditional rectangular plates, the analysis is within 0.5%. 
For simply supported boundary conditions on the parallel ends, the parametric relationships showed that 
increasing the taper of a plate always increased the buckling load when compared to a rectangular plate 
having the same width as the wide end of the tapered plate. There was an interesting interaction between 
boundary conditions and material properties for the plates having clamped ends, however, that caused this 
trend to reverse for certain plate aspect ratios. The clamped end conditions always yielded higher buckling 
loads when compared to plates with the same geometry but with simply supported end conditions. The 
anisotropic parameters were shown to severely decrease the buckling load, and the nondimensional 
parameter p, which includes all of the orthotropic material terms, was shown to have a large effect on the 
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buckling load. Increasing p caused a significant increase in the buckling load, and the decrease in buckling 
load due to anisotropy was shown to be less for higher values of p. The nondimensional parameter a, which 
is the ratio of the longitudinal bending stiffness to transverse bending stiffness, was shown to have a greater 
effect on the buckling load for plates with clamped end conditions when compared to plates with simply 
supported end conditions. For both sets of end conditions there were areas of overlap of the nondimensional 
relationships which indicate areas of flexibility for the designer when choosing stacking sequence and plate 
geometry to achieve a specific buckling load. Contour plots revealed the effects of plate geometry, material 
properties, and boundary conditions on the buckled mode shape. The laminates having their stiff direction 
oriented perpendicular to the loading direction had more half-waves in their buckled mode shape when 
compared to laminates having their stiff direction oriented parallel to the loading direction. The clamped 
boundary conditions were shown to both change the number of half-waves in the buckled mode shape and 
compress the halfwaves toward the center of the plate. The anisotropic parameters were shown to skew the 
buckled mode shapes 

7.3 Conclusions Regarding the Experiments 

A series of experiments were conducted to verify the semi-analytical predictions for a wide range of plate 
geometries and stacking sequences of an AS4/3502 graphite/epoxy composite material. Intrinsic to these 
experiments were test specimens, test fixtures, instrumentation, experimental technique, and determining 
buckling load from the experimental data. The specimen fabrication process, in general, yielded specimens 
having uniform thickness with flat and parallel ends. Specialized test fixtures were required to provide the 
necessary boundary conditions and provide uniform loading across the specimen width without causing any 
out-of-plane bending. The instrumentation was required to verify uniform loading and provide adequate 
information to accurately determine the buckling load. The experimental technique incorporated special 
steps to ensure alignment to provide uniform loading and prevent out-of-plane bending. This combination of 
special test fixtures and experimental technique provided reasonably uniform loading on the specimens. The 
back-to-back strain gages revealed very little out-of-plane bending prior to buckling for all cases, while the 
displacements measured by right and left DCDT’s agreed within 2%, in general, and within 5% for the worst 
cases. The experimental buckling loads were determined by finding the intersection point of straight lines fit 
by hand through the prebuckling and postbuckling portions of the load vs. end-shortening relations. This 
technique worked quite well for the majority of the cases. Some difficulties were encountered, however, 
with the very stiff laminates because there was very little change in slope between the prebuckling and 
postbuckling portions of the relations. 
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The equation describing the prebuckling plate stiffness (slope of the prebuckling portion of the load vs. end 
shortening curve) is slightly more complex for tapered plates than for rectangular plates because the width 
varies with position along the length of the plate. In general, the prebuckling stiffness measured by the 
experiments was within 7% of the analytical prediction. The only exceptions to this were the [0J 24 laminates, 
which had disagreements between 27 % and 35%. With the exception of the [0] 2 4 laminates, therefore, the 
experiments verified that the analysis adequately predicts the prebuckling response of tapered composite 
plates. 

In the baseline comparison with aluminum test specimens, the difference between experimental and 
analytical buckling coefficients were generally within 4%. This verifies that, when material issues are 
excluded from the discussion, the semi-analytic analysis handles the tapered geometry quite well. 
Investigating both [0] 24 and [90] 2 4 laminates as limiting cases for the orthotropic parameters showed 
disagreements less than 8% for the [90] 24 laminates, while the [0] 2 4 laminate had disagreements as large as 
32%. The angle-ply and quasi-isotropic specimens incorporated all of the orthotropic, anisotropic, and 
geometric parameters. In 75% of these cases, the experimental buckling loads were within 13% of the 
analytical predictions. In general there was better agreement between the analysis and the experiments for 
composite plates having their stiff axis oriented perpendicular to the direction of loading. With the exception 
of the [0] 24 case, the experiments verify that the semi-analytic analysis adequately predicts the buckling load 
for tapered composite plates. 

Comparing the shadow moir6 fringe patterns with the calculated buckled mode shapes verified the skewed 
buckled mode shapes due to anisotropy that is predicted by the analysis. However, the semi-analytic analysis 
is not always able to correctly predict the number of halfwaves in the actual buckled mode shapes. 

Finally, the potential sources of error that could possibly contribute to any discrepancy between the observed 
and analytical results were outlined. 

7.4 Recommendations for Future Work 

This research is fairly complete. However, there are logical extensions that can be made. Specifically, a more 
accurate prebuckling force resultant distribution can be incorporated into the analysis. This would require a 
moderate amount of effort. This addition should enhance the results for highly orthotropic plates with large 
tapers that have their stiff direction oriented parallel to the direction of load, and perhaps it would improve 
the prediction of the number of halfwaves in the buckled mode shape. Additionally, this analysis could be 
expanded to include a postbuckling analysis. This postbuckling analysis would require the previously 
mentioned refined prebuckling analysis and would require a significant amount of effort. Finally, this 
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approach could be expanded to incorporate additional geometries, such as plates having a tapered planform 
area with cutouts. These cutouts could be typical, such as circles and rectangles, or could continue in the 
spirit of some lack of symmetry, and be trapezoids, or perhaps even tapered ellipses. 
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Appendix A 

Evaluation of Boundary Integral 


The boundary integral that results from integrating the coefficient of D 66 by parts twice using Green’s 
theorem is 


dA 

The terms n and n are the x and y components of the inplane outward normal vector of the plate boundary. 

x y 

Consider a parametric curve f(S) in a cartesian coordinate system. 



-i- F n 

3xv5yy y dyVdyV 


}ds- 


(A.l) 


Each point on the curve has associated normal and tangent unit vectors n and s where 


fi x = n J fi y = n yJ 

§ x = § x> § y = § yj 


(A.2) 
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The chain rule for differentiation yields the following relations 

df 3f 3s 3f 3n df 3f df df 

— = 1 s + — n = — n n a 

3x 3s3x 3n3x 3s x 3n x 3n x 3s y 

df _ 3f3s ^ 3f 3n _ ft + df _ df df 

3y 3s3y 3n3y 3s y 3n y 3s x 3n y 

df 3f 3x 3f 3y 3f df df df 

— 1 = s + — n = — n n c 

3s 3x3s 3y3s 3x x 3y x dy x dx y 

Using part c of equation (A. 3), the integrand of equation (A.l) can be expressed as 

3w 3 fdw^ 

~3x 3sldy J’ 

Expanding the integrand of (A. 1 ) using parts a and b of (A.3) yields 


< 

_< 


"3w A 3w 
n. n 


.ds y 3n x Jds 3s 


3 r dw 3w 
{r-n,+ — n w } = 


3n 


"3w 3w 

L^ ny “*; n ’ 



3 w 3 2 w 


dn. 


?dwd 2 w 

= n n + n h n — 

x y ds 3 2 y ds dnds y 3s 


^ n x3w ^ n y3wl 

n + + — — } 

y 3s 3s 3s 3nJ 

fdwV ^ n y 3w3w 
V3s ) + "'ds ds d^ 


2 3 w 3 2 w dwd 2 w ^ n x3w3w ^ n y f3wV 

n n n n n — - — . 

^ n 3s 2 y 3n 3nds 3s 3n3s x 3s vdn/ 


Using 


3w3w _ 3 Tdw3w 
3n3n3s 3nLdn3s_ 


d 2 w3w 

dn 2 ^ s 


The integrand can be expressed as 


n n 


f 2 2 

| 3w3 w 3 w3w 3 “3w3w" 
1 3s 3 s 2 3 n 2 3 s 3nL3nds_ 


2 2 
2 3 w 3 w 2 3 w 3 w 

4- n n 

y 3s 3n3s x 3 n 3 s 2 


^ n x 3 w 3 w ^yfdwV fd w> | 2 ^ n y 3 w 3 w 

x 3s 3n3s x 3s \3n J y 3s Ids J y ds 3s 3n 


For constant w everywhere on the boundary, 


(A.3) 


(A.4) 


(A. 5) 


(A. 6) 


(A.7) 
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3w 

— =o 


along the boundary curve. For simply supported or clamped boundary conditions on all edges of the plate, 
therefore, the integrand simplifies to 


3n y pwV 
~* x ds l Bn) 


5W ^ n y 

This will equal zero if — = 0, — = 0, or n x = 0. 
dn os 

Consider a plate of arbitrary geometry having area A and boundary dA . Each point S on the boundary has an 

associated infinitesimal arc length ds with outward normal n (S) which is a function of position and has 

components in the x and y directions n ( S) = n(S)i + n(S)j . As the boundary of the plate is traversed, 

dn y 9w 

the orientation of n changes, so — y * 0. If the plate is clamped on all edges, — = 0, but for simply 

3w . n * . 

supported edges, — * 0. And finally, n^0, except for the two points at which n is aligned with the y axis. 

dn 
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Now consider a polygonal plate with N sides having boundary 9A = 9A, + 9A 2 + . . . + 9A n . The boundary 
integral is now written in a piecewise continuous fashion 




(A. 10) 


Consider a typical edge. Since each edge is a straight line, n x and n y are constant values. This results in 




9s 9s 

Therefore, for polygonal plates having simply supported or clamped edges, 


* 


r 9w 

l 9x 




}ds = 0 


(A-ll) 


(A. 12) 
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Appendix B 
Details of Integration 


The out-of-plane deformation of the plate is expressed as 


where 


N N 

w(^ri) = 5^ 1 wjy, [fj (^, ti) 1 +w°r j [g j (^,ti)]} , 

i= lj = 1 


Vj = cos (fj) , 


= sin (gp , 


and 


(2j-l) 


7CT| 

2b (O’ 


jTtri 

gj b (£) ’ 


b(^) = 1 + e£. 


The derivatives of individual terms of eq. (B.l) are given by 


df; Jtrj 

dt, 2b(£) 3 


a 2 f. 


d? 


1 = 2e 2 (2j - 1) rtT1 


2b (^) 


(B.l) 


(B.2) 


(B.3) 


(B-4) 


(B.5) 


(B.6) 


(B.7) 


(B.8) 
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df, n 

— - = (2j - 1) 

dn 2b(^> 


a 2 f 


— 1 = o, 
dn 2 


a 2 f. 

J_ 


= -e(2j - 1) • 


2b (O' 


= jim 
*“ b(^r 




-> = 2e 2 ™-, 

d$ 2 b(^) 3 

9g j jic 


b($) 


9g j 
2 = °’ 
an 


a 2 g, 


J =- e ^L 


3^1 b(^) z 

Making these substitutions into eqs. (2.88), (2.92), and (2.96) provides 


C‘ E = <D j (^)e(2j-l)-lisin(f ) + O', (^) cos (f ) , 



-<I> j (^)e^YCOs(g j ) +0'^^) sin (gp , 


C 2E = 0"(4) 008(9 +£ (2j- 1 ) ” 

b 





t)sin (9 


2 

-<D i (9e 2 (2j-l) 2 --n 2 cos(9 , 
4b 4 


(B.9) 


(B.10) 


(B.l 1) 


(B.12) 


(B.13) 


(B.14) 


(B.15) 


(B.16) 


(B.l 7) 


(B.18) 


(B.19) 
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cfj° = <D"j (£) sin (g.) - 2e } - 

b 2 




eO. 




.2 2 

-«Dj(^)e 2j — r| 2 sin(g j ) , 
b 4 J 


Njj = -<^(5) (2j - 1) — sin(f ) 
J 2b 


Nj l j° = *,($)^cos( gj ). 


N 2E = (2j - 1) 2 — -^cos (fj) , 

4b 


v 

N 20 = -tf.^J—sintgp, 
b 


M 2E = Oj (^)e(2j - 1) — [sin (f ) + (2j - 1) ^ cos (f)l - (2j - 1) (5) sin (fp , 

J 2o L 2b J iv 


and 


M ?i° = 

b 


cos (gj) + sin (gj) 


+ J ^4>/^cos (gj) 


The total stiffness matrix given by eq. (2.110) is expressed as 


N N N N 


k t (w, 8w) = y y y y uk ee w e +k t ° e w°)8 W e mk t eo w e + k t °°w°)8w°} 

Z-d Lj 1 T ijpq IJ Tjjpq ij' pq T IJpq ij I ljpq ij pq 


i = 1 j = 1 p = lq = 1 


The first term is the series is given by 


K. 


EE 

r ijpq 




0 


(B.20) 

(B.21) 

(B.22) 

(B.23) 

(B.24) 

(B.25) 

(B.26) 
. (B.27) 

(B.28) 
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where (!;) contains the integrals performed across the width of the plate. 

In order to simplify the equations so that the symbolic manipulator MACS YMA can be use most effectively 
for performing the integrations, the equations will be broken down into small parts, the integrations will be 
performed, and then the parts will be reconstructed. The first logical breakdown is to factor the equation into 
the various combinations of the order of differentiation of <|>, namely 

F ee (£) = f 2 V' <D" +f 2 V'.<I>' +f 2 V'.<I> +f' 2 O‘0" +f u 04>’ + f'V<I> 

jq i p jq i p jq » p jq j p jq i p jq » p 

x jq i p jq i p jq i p y 


(B.29) 


Note that the first superscript of fj q represents the order of differentiation of (ft while the second superscript 
represents the order of differentiation of <t> p The fj q are now factored according to the specific integrals over 
TV The Ij q are the integrations over r| across the width of the plate. The first superscript on Ij q signifies which 
of the four possible combinations of trig terms occurs in the integrand according to the following list: l=sin 
sin, 2=sin cos, 3=cos sin, 4=cos cos. The second superscript signifies what power T\ is raised to, i.e. 0=T|°=1, 
1=T|, 2^r| 2 , etc. The fj q are given by 


f 22 = a'l 40 
jq x jq ’ 


(B.30) 


*jq = TT (2q — 1 ) [aeljq -bylj^], 
b 


(B.31) 


f 20 = 
r jq 


~ (2q - 1 ) : 2 [cx 2 e 2 I j 42 - 2aeybI 4 J + b 2 (JlJ] - ^ (2q - 1) [aeljj - byl jj] 
4b b 


(B.32) 


f' 2 = 5S5 (2j - 1) lael’ 1 - byl’?] = f 2 '. 

jq 2 J qj 1 qj qj 

b 


(B.33) 


f j<! = (2j_1) (2q_1) 


(B.34) 


f 10 _ 

jq “ 


(2j - 1 ) (2q - 1) : 2 [aVl* - 3aVb Y I 22 + aPeb 2 ^ 1 - b 3 SlJ°] 

4ab 

(2j- 1) (2q — 1) [ cxeljq — 2bylj” ] , 


(B.35) 
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f? q 2 = -JL ( 2j-i) 2 [aVlJ-2aeb7i; q , + b 2 piJ] [aelj|-l<] = . (B.36) 

4 b 


Jq .3 
b 


% = (2q - 1 ) (2j - 1 ) : 2 [aVl* - 3aVb Y l“ + «Peb 2 I 3 ‘ - b 3 5l 3 q °] 


Jq ^ *6 
4ab 


(2j - 1) (2q - 1) [aelj 2 - 2byI* Q '] = , 


(B.37) 


and finally 


*jq = 1 ~~T^ (2j - 1) (2q _ 1) [ a£ Ijq - 2bYlj q ] 
b 

7i 3 e 


+ JLL (2j - 1) (2q - 1) { a 3 e 3 [ (2q - 1 ) I 23 + (2j - 1) I 23 ] - 3a 2 e 2 b Y [ (2q - 1) 1“ + (2j - 1 ) I“] 
4ab 7 Jq (B.38) 


,22 


22 , 


+ apeb 2 [ (2q - 1)1^ + (2j - 1) I 3 q '] - b J 6 [ (2q - 1)1™ + (2j - 1) I jq ] } 

+ - 2 — ( 2q - 1 ) 2 (2j - 1 ) 2 [ aV I 44 - 4aV b Y I 43 + 2aVb 2 pi 42 - 4<xe8b 3 I 4 ; + b 4 I^] 


? 8 

16ab 


Because the specific integrals over T| are now separated out, they can be easily recognized and computed 
using the symbolic manipulator MACSYMA. The results are given as 


I 40 

jq 


* -b (,k) 


COS 


( 2 j - 1 ) 


7CT| 

2b^)J 


cos 


(2q — 1) 


JlTl 

2b (X) 


Jdri = 0;j*q 

= i ; j = q 1 


(B.39) 


A-r " 

" -b (£,) 


COS 


(2j - 1) 


7TT} 

2b (%) 


•J sin |^(2q - 


1) 


KT\ 

2b(f)J 


dt| 


b 2 ( 2 j — 1 ) (-l) J+q 

Jt(q-j) (q+j- 1) 

b 2 

7t f 2i — 1) ’ J ~ q ’ 


;j*q 


(B.40) 
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=1 


-*(£) 2 

T] COS 


(2j— 1) 


7tT) 

2b(£). 


b($) 

_ 2b 3 (2j - 1) (2q - 1 ) (~l) J+q 
ft 2 (q-j) 2 (q+j- 1) 2 

_ b 3 [7t 2 (2j — 1 ) 2 — 6] 

37t 2 (2j-l) 2 ’ J q ' 


COS 


( 2q — 1 ) — — — 1 dr) 

. H 2b(^)J 


;j*q 


i 12 

jq 


3 r 7cri n r 

= J q sin ( 2j - 1) 2b(^) jsin[(2q- 0 


-b(4) 

2b 3 [ 2 j (j — 1 ) + 2q (q - 1) + 1] (-1) 


7CT| 

2b (%). 


j+q 


ft 2 (q-j) 2 (q+j- 1) 2 


-;j*q 


b 3 [n 2 (2j - 1) 2 + 6] . _ 

>J — Q - 


3ft 2 (2j - 1) 2 


= j n 


ri 3 sin 

b ((,) 


(2j - 1) 


jcq 

2b (^). 


cos 


[(2q — 


1) J-ldtl 

2b(^)J 


b 4 (2q - 1) {q (q — 1) [ft"(q (q- l) -2JG - 1) )~6]} (-l J + q ) 


re 3 (q — j) 3 (q + j — i) 3 

b 4 (2q — 1) {j(j- 1) [ft 2 j (j ~ 1) ~ 18] -6} (-l Hq ). : 


ft 3 (q - j) 3 (q -b- j - 1) 3 

b 4 [ft 2 (2j — 1 ) 2 — 6] . 

= 7 ^ — ;j = q’ 

ft (2j-l) 3 


j*q 


4 

jq 


I 


T) COS 
-b<$) 


(2j - 1) 


7tT) 

2b(^) 


j cos j^(2q - 


1) 


ftt] 

2b(^)J 


dq 


4b 3 (2j - 1) (2q - 1) {q (q - 1) [ft"(q (q-1) - 2j ( j - I ) ) -12]} (-l , + q ) 


ft 4 (q-j) 4 (q+j- 1) 4 
2 - , 


+ 4b 3 ( 2 j — 1 ) (2q - 1) {j (j - 1) [ft‘j(j-l) -12]-6} (~l j4q ) 


ft 4 (q-j) 4 (q+j- !) 4 


j *q 


b 5 { 7r 2 ( 2j — 1 ) 2 (ft 2 (2j - l) 2 -20) + 120 } . _ 

5^(2j-l)‘ ' , “'' 


I 30 

jq 


I 


■*<$) 
cos 


( 2 j — 1 ) 


-b($) 


icq 

2b(^)J 


sin 


(2q — 1 ) 


Jtq 

2b(^)J 


dq = 0, 


(B.41) 


(B.42) 


(B.43) 


(B.44) 


(BAS) 
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41 f rt(4) 

i jq - J ncos 


(2j-D 


-b<5) 


JtT) 

2b O) 


■ J cos ^ 


(2q — 1) 


nil 

2b (^) 


\dr\ = 0, 


I." = f r|sin (2j - 1) — H_1 s in f(2q - 1) — ^Idri = 0, 

Jq J 1 L J 2b(^)J L V 2b(^)J 

-b<5) 


j* r| 2 sin ( 2 j — 1 ) WT ) 1 cos [~ (2q - 1) dr] = 0, 

J L 2b (4) J L H 2b(^)J 


22 . 

I 22 = I i, 2 ™ 

-b(5) 


43 f*®J 

I -l = n cos 

-b($) 


( 2 j - !) cos [( 2 q - 1) WT1 


2b (§). 


2b (^) 


]dt| = 0. 


(B.46) 


(B.47) 


(B.48) 


(B.49) 


Combining the evaluated integrals of eqs. (B.39) through (B.49) with the appropriate coefficients in eqs. 
(B.30) through (B.38) and substituting into eq. (B. 29) yields 


.4>- p d5 + cJ| + 

0 0 0 

+ c 21 fi>' <D" d^ + c’ 1 d^+c^f-dJ'.^ + c” f I®.®" d£ 

qj J > p ^ jqJ b p ^ jq J b 2 p b p 

o o o° o 

+ c‘°f— OO' d^ + c°°f 

qjJ , 2 1 p ^ jq J b 3 P * 


(B.50) 


where 


22 

C Jq 


= o. 


(B.51) 


C 21 = q 2 £(2j - 1) (2q — 1) (~l) i+q 

Cjq (q-j)(q+j-i) 


(B.52) 


.20 = « 2 e 2 ( 2 j — 1 ) (2q - 1) (— l) j+q [2q (q — 1) +2j(j-l) + 1] 
jq 2(q— j) 2 (q+j — l) 2 


(B.53) 
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ii = 2aV(2j-l) (2q - 1) (-1) J + q [2g(g-l) +2j(j-l) +1] 

Cjq (q-j) 2 (q+j-D 2 

c? q ° = -aV (2j — 1) (2q— 1) (-1)- j + q {q 2 (q - l) 2 [n 2 ( (2q - 1) 2 - 8j (j - 1 ) ) -8] 

+ 2j(j — l)q(q-l) [7t 2 (2j 2 — 2j - 1) -36] 

— 22q (q - 1) +j(j-l) [ (it -4) (Jt + 4)j (j - 1) -26 ] +-6 } / (4 (q - j) 3 (q +j - 1) 3 ) 
-K 2 pE ( 2 j — 1 ) (2q- 1) 3 (-1 ) j + q / (4 (q -j) (q+j-1)) - 


(B.54) 


(B.55) 


c oo _ it 2 pe 2 (2j-l) (2q - 1) (-l) i+q [2q(q- 1) + 2j(j-l) + 1] 2 
Cjq 4(q-j) 2 (q+j-l) 2 

+ aV(2j-l) (2q - 1) (-l) j + q [2q (q- 1) + 2j(j- 1) + 1] (q(q- 1) 

[q (q — 1) [K 2 (2q(q-1) -2j (j - 1) +1 ] -4] 

-2j(j-U [n 2 (j 2 -j+l) +44] - 24] + j(j-l) [ j ( j — 1 ) [tc 2 (2j 2 - 2j + 1) 
— 4] —24] - 6}/(4(q-j) 4 (q+j— l) 4 )- 

For the cases when j equals q then the c qq terms become 


22 2 
= a , 

qq 

0.57) 

21 -.2 
= a e, 

qq 

(B.58) 

" a£ [7i 2 (2q-l) 2 + 6] - n 4 P (2q-l) 2 , 

(B.59) 

2 2 


-^ln 2 (2q-l) 2 + 6) , 

(B.60) 

" aE [7n 2 (2q — 1) 2 + 6] — (2q — 1) 2 , 

(B.61) 


and 
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2 4 

c 00 = — {7i 2 (2q- O 2 [37t 2 (2q - 1) 2 + 140] + 120 } + 
qq 240 

+ -^L(2q-l) 4 . 

80a 2 


^i^(2q-l) 2 [7C 2 (2q-l) 2 + 6] 

24 


(B.62) 


The integrands in equation (B.50) are functions only of J; and will be different for the clamped or simply 
supported end boundary conditions and will be considered later in this appendix. The term can be 

expanded, the integrations performed, and terms recombined in the same manner as the K T ..^ term ’ 
specifically 


K, 




OO | "cOO 

r ijpq 


where 


(B.63) 


F ijpq (0 = hj q V p + h>- p p + hj>" + h’ 2 ^'^” p + h^vo-p + h^'^p 




(B.64) 


The hjq are given by 


s ■ «*>;■ 


(B.65) 




(B.66) 


u 10 - 
h jq “ 


• <B ' 67) 


l.;, ! = -^[ael”-b T j”l -I.” 
b 


(B.68) 


h » = 4j^m [aE] «-2 byj£l. 


(B.69) 
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and 


h“ = la y ,» _ aaVbyjJJ + ap E b J jJ' - b J 5j”] 


jq .6 
ab 


jq 


_47 L 0^q [a£j 42_ 2byJ 4,] ) 

b 5 


jq J 


h ” - +^,«< , -brl”) - <■ 

b b 


. }^h [a’^-aaVbrJ^aM^-b’Sj”] 
ab 

_2 2 . 


_4na£jq « _ 2b j4. j = h .o 

5 L jq * jq J <U * 


-OO 4rc 2 ae 3 jq . x 42 „ x 4i 


1 uu 

h jq = 


. t a£J iq- 2b ^] 


{aV [qJ^+jjJ 3 ] -3aVbytqJ 3 q 2 +jjJ] 

ab 

+ aPeb 2 [qJ 2 J , + jJ 2 ’]-b 3 5[qjJ + jj“] } 

+ [ aYf 4 - 4a 3 e 3 b Y j; 3 + 2aV b 2 pj; 2 - 4ae5b 3 J j 1 q 1 + b 4 jJJ]. 

a z b 


The Jj q are the integrations over r\ across the width of the plate, given as 


b(t) 


= 1 


sin 


j7UT_ 

b(^) 


sin 


qrni 

Lb (5) 


-b<$) 


J dT| = 0;j*q 
= b ; j = q 


J 2, = 

jq 


b($) 

J 

-b(£) 


rjsin 


j tzt \ 1 f qTcri 

Lb(^)J Lb (5). 


dri = 


2b 2 j(-1) i + q 

it(j-q) (j + q) 


; j *q 


2qJt 


;j = q 


(B.70) 


(B.71) 


(B.72) 


(B.73) 


(B.74) 


(B.75) 
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b (£) 


1 


r| 2 sin 


' jrcri 

b(^). 


sin 


qnri 


8b 3 jq(-l) J+q . 

drt = . JH \ - !J * q 


-b<5) 


b (^)J Jt 2 (q-j) 2 (q + j)' 


b 3 (27t 2 q 2 — 3) 

6n\ 2 ;J_q 


J 42 

jq 


b(5) 

-bU) 


COS 


jwn 

Lb(^) J 


cos 


q*n 

■b(^) 


4b 3 (q + j) 2 (-l) j+q 

dT| = ~1 2 2 ’ J q 

^ 2 (q-j) 2 (q + j) 

b 3 (27C 2 q 2 + 3) . 

= — 2 ;J=q ’ 

6nq 2 


b (£) 


J 


3 

T) COS 


-b(^) 




sin 


qrcr) 

b($). 


dt) 


2b 4 q{q 2 [7C 2 (q 2 -2j 2 ) -6] + j 2 (tt 2 j 2 - 18) } (-1) j +q 
Jt 3 (q-j) 3 (q+j) 3 

b 4 (27t 2 q 2 -3) . 

= Z7^ ;J = q ’ 

4Jt q 


;j*q 


b(5) 


j ' 4 

jq 


f 4 . 

‘ jnq 

sin 

■qjiq“ 

I T| sin 

Lb(^) J 

Lb(^)J 


dq 


-b(^) 


16b 5 jq { q 2 [ 7t 2 (q 2 — 2j 2 ) - 12] +j^Jcy-12)H :: l) 


.2 . 2.2 


,j + q 


Jt 4 (q-j) 4 (q+j) 4 

b 5 [27t 2 q 2 (7t 2 q 2 -5) + 15] . = 

10n 4 q 4 


;j*q 


b()j) 


j 


20 


jq 


r . r ■ 

= Sin — — 

J Lb(^) J 




Lb(^) J Lb(^)J 


dq = 0, 


-b(^) 

b(4) 

Jj q ' = J -n sin 

-b<$) 

b<£) 

I 

-b($) 




.b(?)J Lb© 


dq = 0, 


f „ 

r jitq " 


r q*n“ 

1 T| COS 

Lb(^)_ 

COS 

Lb(^)J 


dq = 0, 


(B.76) 


(B.77) 


(B.78) 


(B.79) 


(B.80) 


(B.81) 


(B.82) 
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and 


b (£,) 




.32 | 2 f j^T| 

J jq = I t| COS | 


Lb(^) 


sin 


b(^)_ 


dri = 0, 


-b(^) 


J ? 3 = 

jq 


b (4) 

J 

-b(4) 


r| 3 sin 


1 sin f I” 1 


Lb(yJ Lt>($>. 


dr| = 0. 


(B.83) 


(B.84) 


Again, the first superscript on Jj q signifies which of the four possible combinations of trig terms occurs in the 
integrand according to the following list: l=sin sin, 2=sin cos, 3=cos sin, 4=cos cos. The second superscript 
signifies what power T| is raised to, i.e. 0=r}°=l, l=r|, 2=r| 2 , etc. 

Combining eqs. (B.65) through (B.84) and substituting into eq. (B. 64) yields 


k c «> ■ "Sji*- i®" . v? 


+ 


d 2 ' iV*' d^ + d" f -O'.O' d^+d'T-O'.O d^ + d^f-M" d^ 

<Uj » P ~ J<3j ^ 1 P ^ JQJ ^2 1 P ^ <Uj ^ 1 P ^ 


(B.85) 


+ d'° f— ' d^ + d°° f — 0,0 dE, 

<U J ,2 1 P ^ ,3 1 P ^ » 


0 


0 


where 


< - «. 


(B.86) 


d 2i = 4a 2 £jq(-l) J ' t ' q 

Jq (q -j) (q+j) 


(B.87) 


d 20 _ 4aVjq(q 2 +j 2 ) (~l) J+q 
Jq (q - j) 2 (q + j) 2 

d n _ 16a 2 £ 2 jq (q 2 + j 2 ) (~l) J + q 

Jq (q-j) 2 (q+j) 2 


(B.88) 


(B.89) 
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and 


d io = 47t 2 pejg 3 (-l) ) + q 

jq (q-j) (q+j) 

4a 2 e 3 jq { Tt 2 (q 6 - 2j 2 q 4 + j 4 q 2 ) + 4j 4 - 1 8j 2 q 2 - 10q 4 } (-1 ) 3+q 

(q -j) 3 (q+j) 3 


,00 = 47C 2 |3e 2 jq(q 2 +j 2 ) 2 (-l) 

,s2,_ . ; ,2 


j + q 


+ 


(q-j) (q+j)‘ 

4a 2 £ 4 jq (q 2 + j 2 ) { ( q 6 — j 2 q 4 — j 4 q 2 + j 6 ) - 2j 4 - 44j 2 q 2 - 2q 4 } (~1) J q 

(q — j) 4 (q + j) 4 


If j equals q then the d qq terms become 


-e - “ , »- 


- <*. 


2 2 


d 20 q = -^(27l 2 + 3)-n 2 Pq 2 , 


d 


11 

qq 


2a 2 e 2 

3 


(2rc 2 q 2 + 3) , 


d'° = 
qq 


a 2 £ 3 


(27i 2 q 2 + 21) -PetcV. 


2 2 


and finally 


d“ = UlL [27t 2 q 2 (3n 2 q 2 + 35) +15] + ^-^q 2 (2jt 2 q 2 + 3) + ^-. 


The K t oe term can be expanded, the integrations performed, and terms recombined in the same 

T ijpq 

namely 


K 


°E _ 
T ijpq 



(5) dE, 


(B.90) 

(B.91) 

(B.92) 

(B.93) 

(B.94) 

(B.95) 

(B.96) 

(B.97) 

manner, 

(B.98) 
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where 


F oe (?) = r 22 <t" <D" + r 21 <&" 4>' + r 2 V' <D + r 1 V +r 1, 4> a>' H-rV.O 

npq v ^ ; jq i P JQ » P JQ i P JQ i P jq > P jq » P 


+ r 02 O O" + r 01 O<I>’ + r°V<D 
jq i p jq > p jq * p 


The rjq are given by 


22 2 t 20 
F jq = ° L W ’ 


r£ = (2q - 1) [otel-jq “byLjq] , 

b 


20 

r jq ~ 


( 2q - 1 ) 2 [ a 2 e 2 L 22 - 2«E Y bL 2 ; + b 2 pL 2 J] 
4b 


(2q - 1 ) [ a£ Ljq _ byL. 1 ®] , 
b 


12 2jraj _ _ 4 i . T 40i 

, « = -rr |atL ir b » L ). 1 ' 


11 271 CX£. . v r 32 x 31i 

r j q = — J (2q - 1) [aeL jq - 2byL J , 


jq J 


10 11 i \ 2 r 3 3y 4 3 o 2 2. T 42 o . 2 T 41 .3^x40, 

r ; _ = -j (2q — 1) [a e L jq -3a £ b Y L jQ +<*P£b L ja -b SL ia ] 


Jq 2 ab 6J 


jq 


jq 


jq J 


27 t 2 a £ 2 


5 j(2q-l) [a£L 22 -2b Y L 2 ( |] 


r ° 2 = -^L [a 2 e 2 L 22 - 2a£ Y bL 21 + b 2 pL 2 °] + ^!j [ck-L* 1 - byL^] 
jq 4 L jq * jq K jq J .3 J jq jq J 

b b 


r?. 1 = — — -r j 2 ( 2q - 1) [aVl^ 3 - BocVbyLJ 2 + aPEb 2 ^ - b 3 8L? 0 °] 
ab 


jq 


jq 


jq 


jq J 


27t 2 ae 2 


5 j(2q-l) [a£L 32 -2b Y L 3 |] , 


(B.99) 


(B.100) 


(B.101) 


(B.102) 


(B.103) 


(B.104) 


(B.105) 


(B.106) 


(B.107) 
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and 


r oo = . (2q _ , } f aeL £ _ 2 byL 3 ’ ) 


jq 


jq J 


+ —J (2q - 1) {a¥ [2jL“ - (2q - 1) L 43 ] - 3aVby[2jL j ' q 2 - (2q - 1) L 42 ] 


2ab 7 


+ apeb 2 [2jLj q - (2q - 1) Lj] - b J 5 [2jL‘ q u - (2q - 1) L£] } 

^-j 2 (2q - 1 ) ' 2 [ oc 4 e 4 L 24 - 4a 3 E 3 byL 23 + 2a 2 e 2 b 2 pL“ - 4ae5b\f q + b 4 L 2 q °] . 


jq J 


jq 


a„ 2 h 8 


The Lj q are the integrations over r| across the width of the plate, given as 


btf) 


L 10 = 
jq 


r • 

‘ jTCTl " 

i . r 

1 sin 

J 

Lb(^) J 

i 

c 

7 a 


(2q - 1) 


7tri 

2b ($) 


Jdri = - 


8bj (— 1) 


j+q 


lt(2q-2j-l) (2q + 2j — 1) 


-b(^) 


b(5) 


L 

jq 


■ J 


T|sin 


' j*n i ; 
■b(5)_ 


sin f (2q - 1 ) WT t -~| dt] = 0 

L 4 2b(t)J 


-b(4) 


jq 


b(^) 


-b(4) 


f « 2 

jJCTl "1 

r 

I T) sin 

Lb (!;) _ 

sin 


(2q — 1) 


ItTl 

2b (5). 


dr) = 


-8b 3 j { 8q (q — 1 ) [it 2 (2q 2 - 2q - 4j 2 + 1) - 12] + 8j 2 [li 2 (2j 2 - 1) -4] + it 2 -24} (- 

Jt 3 (2q — 2j - 1) 3 (2q + 2j — 1) 3 


b($) 


f t| 3 sin s i n (2q-l) ^ dq = 0, 

J Lb($)J L' 2 b(^)J 


Ljq = | n 3 sin 


-b($) 



b(b 



II 

8 . 2 ; 

J 

1 sin 

1 1 

CT 

v/TC j 

1 1 

cos 


(2q — 1) 


ItT) 

2b (0 


•Jdri = 0, 


-b (4) 


b(0 


L 21 = 
jq 


f 

j 7tr| 

‘ 

1 rjsin 

Lb(^)_ 

cos 


(2q — 1 ) 


7tri 

2b (0. 


32b 2 j (2q — 1) (-1) 


dr| = 


j + q 


-b(£) 


it 2 (2q — 2j— l) 2 (2q + 2j- l) 2 


(B.108) 


(B.109) 


(B.110) 


(Bill) 
l)j + q 


(B.112) 


(B. 113) 


(B.114) 
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L 22 

jq 


b(^) 



r 2 • i 

J ^ sin 

jTITl 


Lb (£) J 

vvO 


( 2 q - 1 ) 


JtT| 

2b (0. 


dri = 0, 


-b ((,) 


b<0 


■J 


f -n 3 * 

' j*n " 

r 

I T| sin 

Lb(^)J 

COS 


(2q — I) — ldri = 

2btf)J 


-b(S) 


96b 4 j(2q- 1) { 8q (q — 1 ) [k 2 (2q 2 - 2q - 4j 2 + 1) -8] + 8j 2 [rt 2 (2j 2 - 1) -8] +TC 2 - 16} (- 

7l 4 (2q-2j-l) 4 (2q + 2j-l) 4 


b (£) 


L 24 * f n 3 sin [jH-l cos [ (2q - 1) -^-Idri = 0, 

jq J Lb($)J L V 2 b(^)J 


-b<$) 


b(5) 


L-' = I 11 cos 


jq 


1 


jrcil 

Lb(^) J 


sin 


( 2 q - 1 ) — — ^ — ~| dr| 

2 b(^)J 


-b($) 


8b 2 [(2q— l) 2 + 4j 2 ] (-l) j+q 
it 2 (2q — 2 j — 1 ) 2 (2q + 2j — 1) 2 > 


b(5) 

L 32 = f t| 2 cos ~j sin f(2q — 1) Idr) = 0 

Jq J Lb(^) J L v 4 2 b(^)J 

-b(4) 


b(^) 


40 

"jq 


f 

1 

-=$ 

i r 

J cos 

Lb(^) J 

o 

o 

on 

1 


( 2 q - 1 ) 


jcq 

2b 


■]d n = - 


4b (2q — 1) (~l) j + q 
71 (2q — 2 j — 1 ) (2q2j — 1) 


-b(4) 


b($) 


L 4 ' = I T 1 COS 


jq 


I 

-bU) 


■jinn 

-b (^)J ' 


cos 


[ (2q - ,, irik>’’ = 0 ’ 


b($) 


42 

J jq 


f 2 

j«n 


1 rj cos 

Lb(^) J 

COS 


(2q — 1 ) 


7tt| 

2b(^). 


dri = 


-b(S) 


4b 3 (2q - 1) { 8q (q — 1) [Jt 2 (2q 2 - 2q - 4j 2 + 1) -4] + 8j 2 [it 2 (2j 2 - 1) -12] +K 2 -8} (- 

n 3 (2q - 2j - 1) 3 (2q + 2j — 1) 3 


(B.115) 


(B.116) 

l) i+q 


(B. 117) 


(B.118) 


(B.119) 


(B.120) 


(B.121) 


(B.122) 

l) j + q 
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and 


b(S) 


43 = f r| 3 cos cos (2q-l) 

jq J Lb(^) J L v 


rcr| 

2b(^j. 


dii = 0. 


-b(4) 


Combining eqs. (B.65) through (B.84) and substituting into eq. (B. 64) yields 


K C ® = j b®" ,®" ,<« + 4' J®' .®V« - <ffj‘ j®" , Vt 

0 0 0 

+ *■ ,<* + «,"j" J®', v* *cf J®.® - 

oooo 

+ e°‘ f— <&=<!*’ d^ + e“ f — «D d> 
j qJ ^2 1 p ^ jqJ ^3 1 p ^ * 


where 


e 2i = 8ayj (2q — 1 ) (-l) j+q = _ e i 2 

ejq (2q — 2 j — 1 ) (2q + 2j— 1) jq ’ 

e 20 = 8ayej (2q - 1) [ (2q - l) 2 + 4j 2 ] (~l) J + q = ^ 

Cjq (2q — 2j-l) 2 (2q + 2j-l) 2 ,q ’ 

e ii = 32ayej(2q-l) [ (2q- 1 ) 2 + 4j 2 ] (-l) j+q , 
jq (2q — 2j — 1) 2 (2q + 2j - l) 2 

jo 27t 2 8j(2q-l) 3 (-l) ,+q 
C j q a(2q — 2j — 1) (2q + 2j — 1) 

+ 2aye 2 j (2q — 1) (4q(q-l) (4q(q-l) [3tt 2 (4q 2 - 4q - 8j 2 + 3) -8] 

+ 48j 2 [Jt 2 (j - 1) (j+1) -6] + 9tt 2 - 16 ) 

+ 8j 2 f2(3jt 2 - 16) j 2 - 3tt 2 - 361 +3 ji 2 -8 } (-1) j+q / f (2q - 2j - 1) 3 (2q + 2i - 1 ) 


(B.123) 


(B.124) 


(B.125) 


(B.126) 


(B.127) 


(B.128) 


(B.129) 
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8 JT 2 5j 3 (2q — 1) (— i) j+q 
a(2q-2j-l) (2q + 2j— 1) 

2 2 2 2 2 (B.130) 

-8ayej (2q - 1) {8q (8q - 1) (2q (q - 1) [3% 2 j 2 - 4] - 3j 2 [n 2 (2j - 1 ) (2j + l) +12] -4) 

+ j 2 (8j 2 [37t 2 (2j 2 -l) -4] +3tc 2 -72) -4 } (-1) j + q / [ (2q -2j - 1) 3 (2q + 2j - 1) 3 ] > 


c oo_ 27C 2 5ej 3 (2q - 1) (4q 2 -4q + 4j 2 + l) 2 (-l) j + q 
Jq a(2q-2j-l) 2 (2q + 2j-l) 2 

-2aye 3 j (2q - 1) ( (2q - l) 2 + 4j 2 ) {4q(q-l) (4q(q-l) J3J) 

[37C 2 (4q 2 -4q+4j 2 + 3) - 8] -8j 2 [3n 2 2j 2 + 1+88] +9 tc 2 — 16) 

+ 4j 2 (4j 2 [ 3 7t 2 ( 2j - 1) ( 2 j + 1 ) -8] - 3it 2 - 176) + 3 jc 2 - 8 } (-l) j+q 
/ [ (2q - 2j — 1) 4 (2q + 2 j - 1) 4 ] 

Because the integrals associated with these terms have no degenerate cases for j equals q, special 
consideration need not be given for e qq terms. 

Due to the symmetry of the total stiffness matrix the following relationship must be true 

K t eo (§) = K t oe . (5) . (B.132) 

1 MPq ^pqu 


This implies that 


where 



(B.133) 


F e0 (£) = r 22 0" O" + r 12 0" O' + r 02 O" O + r 21 0’ O'' + r n O' O' + r 0I O'O 
upq^ 7 qj »pqj * p qj ip<u»pqj«pqj‘p 


+ r 20 OO" + r ,0 O O' + r°°0 0. 

qj i p qj i p qj » P 


(B.134) 


At this point the integrations over T| are complete and the integrations over £ must now be performed. The 
integrands of these integrals have vastly different functional forms depending on the boundary conditions on 
the ends of the plate. For simply supported boundary conditions on the ends of the plate, and <t> p are given 
by 
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4>j = sin O 1 ^) 
<l>p = sin (pJt£) 


(B.135) 


Making this substitution, the integrands appearing in eqs. (B.50), (B.85), and (B.124) are given by 


1 


-hM = eM 11 + M 10 

Vi Y p ^ ip 

b 3 


>p 


(B.136) 


f = 2e (eM?p + M™) + Tip (£ 2 M 2 p 2 + 2eM 2 p ' + M J) , 

J o b 

J = - 7t 2 P 2 (e 3 M? 3 + 3e 2 M i 12 + 3EM” + M‘ p °) + 2e 2 (eM? p ' + Mj p °) 
0 

+ 47te P (eXp + 2eM fp + <) 


(B.137) 


(B.138) 


J', 

o 


= 2e (eM jp + M^) + ni (£ 2 M 2 p 2 + 2EM 2 , 1 + Mj) 
b 


Ji 

o 


r 


^i^pd^ = 4 e 2 (EM' p + M- p ) 


+ 27t£ [p (E 2 M 2 p 2 + 2EM 2 ; + M 2 p °) + i (£ 2 M 3 p 2 + 2EM 3 ; + M 3 °) ] 

+ n 2 ip(£ 3 Mj + 3e2 < + 3eM X + M ? ’ 


(B.139) 


(B.140) 



o 


- 2k 2 e P 2 (e 3 M; 3 + 3e 2 m’ 2 + 3EM' p + m; p °) + 4 e 3 (eM” + M? p °) 


+ 8Jt£ 2 p (e 2 M 22 + 2EM 2 p +M 2 p ) 


(B.141 ) 


7t 3 ip 2 (£ 4 M 3 p + 4 e 3 M 33 + 6e 2 M 32 + 4eM 3 ; + Mj) + 2jt£ 2 i (E 2 M 3 p 2 + 2 E M 3 p + M 3 p °) 


+ 4n z £ip (e j M 43 + 3e‘M£ + 3EM- + M~) 


2, .42 


.41 


.40, 


= - * 2 i 2 (e 3 Mj p 3 + 3£2m !p + 3eM 'p + M 'p 0) + 2e ' (eM ip + M ip° } 

0 

+ 4it£i (£ 2 M 22 + 2EM 2 p + M 2 p ) , 


(B.142) 
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f'V<M5 = - 2n 2 ei 2 (e 3 M. 13 + 3e 2 m! 2 + 3eM'' + m'°) + 4e 3 (eM'' +M.'°) 

1 |^I Y p ~ x Ip Ip Ip ip 7 v ip ip 7 

0 

+ 8tt£ 2 i (e 2 M | 2 p + 2£M 2 p + M 2 p ) 

- TtVp (E 4 M 3 p 4 + 4 e 3 M 33 + 6e 2 M 3 p + 4eM 3 ; + M 3 p °) + 2ne 2 p (e 2 M 3 p 2 + 2 eM 3 
+ 4n 2 eip (e 3 M 43 + 3e 2 M 42 + 3eM 4 p ' + M J) • 


(Xte = Jt 4 i 2 p 2 (e 5 M* 5 + 5 e 4 m! 4 + 10e 3 m' 3 + 10e 2 m! 2 + 5eM'' + M*°) 

I jyi'P ” r v ip ip ip ip ip ip/ 

0 

-2eV (i 2 + p 2 ) (E 3 M i 13 + 3E 2 M‘ p 2 + 3 eM* p + M? p °) 

-47C 3 Eip[E 4 (iM 2 J + pM 34 ) +4E 3 E 4 (iM“ + pM 33 ) + 6£ 2 (iM“ + P M 3 p 2 ) 
+ 4E(iM 2 p +pM 3 p ) + (iMj + pMj) ] 

+ 8 tie 3 [E 2 (pM 22 + iM 32 ) + 2 e (pM? p ‘ + iM 3 p ) + (ipM* + iM 3 p °) ] 

+ 16K 2 E 2 ip (e 3 M 43 + 3e 2 M 42 + 3eM 4 ; + Mj) . 


The M ip are the integrals over E, as given by 


] 

m'° = f sin (i 7 t 2 ;) sin (ptt^) d£ = 0 ;i *p 
o 

1 . 

= 2 ;,=P * 


Mj" = f^sin(iTt^) sin(p 7 t^)d^= -^;i*p 

J Jt 2 (p-i) 2 (p + i ) 2 


~ 4’ 1 = P 


M ip = |V sin sin (P^) = 7 4ip ^ 54 P 

J 7 t 2 (p-i) 2 (p + i ) 2 


2n 2 p 2 - 3 
127t 2 p 2 


;< = p , 


(B.143) 
+ M 30 ) 

ip y 


(B.144) 


(B.145) 


(B.146) 


(B.147) 
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m ' p 3 = J*£, 3 sin (iTt^) sin(prc^)d^ 


6ip [p 2 ( (tc 2 (p 2 - 2i 2 ) -4) (— 1) l + p + 4) + i 2 ((ni-2) (7ci + 2) (~l)" p + 4)] (B.148) 

rc 4 (p-i) 4 (p + i) 4 


2 2 ~ 

71 p - 3 . 

= — ;i = p , 

o 2 2 r 1 

8tc p 


m' 4 = jYsin ( i7C ^) sin (P^) d £ 


= 8ip [p 2 (TC 2 (p 2 - 2i 2 ) - 12) +i z (rcV-I2)] (-l) IT|, . j 


.2 . 2.2 


fc i + P 


(B.149) 


rt 4 (p-i) 4 (p + i) 4 

2kV(**P 2 -5) +15 . 

^ ;i = P * 

20jc 4 p 4 




= £ sin (irc£) sin (prt£) 


= 10ip { [p 4 (Jt 4 p 4 - 24n 2 p 2 + 72) + 4i 2 p 2 (-n 4 p 4 + 6ji 2 p 2 + 60) +Ji 2 iV(67t 2 p 2 -4jc 2 i 2 + 2^ 15Q) 


+ i 4 (7t 4 i 4 -247c 2 i 2 + 72) ] (-l) i+p - 12(6p 2 + 2i 2 ) (p 2 + 3i 2 ) } / (n 6 (p -i) ° (p + i) °) ;i * p 


6\ 


2 2 2 2 - 15 . 1 A c 

tc p (271 p ) +45 . 

— — ;i = p 

24rc 4 p 4 


M 


l 

r if r-n l+p - n 

20 - I (sin (iJC^) cos (pit^) d£) = — 7 — — ;i*P 


■p 


o 


it(p-i) (p + i) 
= 0;i = p • 


(B.151) 


I 

M ip = f^(sin(iJi^)cos(pn^)d^) = 

,p J jt(p-i)(p + i) 


i(-D 


l+p 


; 0*p) 


(B.152) 


4rcp 


;i = p , 
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Mf p 2 = JV(sin cos (pn^)d^) 
o 

= i[p 2 ((?t 2 (p 2 -2i 2 ) -6) (~l) i4p + 6) + i 2 ((rc 2 i 2 -2) (-l) i + p -f2) ] . 

7T 3 (p — i) 3 (p -4- i) 3 

1 . 

= -— ; | = p . 

4tcd 


Vlfp = J ( sin (i nfy cos (put,) d£) 

0 

i [p 2 (ti 2 (p 2 - 2i 2 ) - 18) - 18 + i 2 (n 2 i 2 
n 3 (p-i) 3 (p + i) 3 


- 2Jt 2 p 2 - 3 . 

= ;i = p 

8jc 3 p 3 


6) ] (-l) i + p 


;i/p 


M 24 = J ( sin cos (P 7 ^) d £) 

0 

= i{ [p 4 (TC 4 p 4 -36n 2 p 2 +120) + 4i 2 p 2 (- n 4 p 4 + 15ic 2 p 2 + 60) + 2n 2 i 4 p 2 (3ic 2 p 2 -2itV • 

+ i 4 (7t 4 i 4 - 127c 2 i 2 + 24) ] (-1) 1 +p - 24 (5p 4 + 10i 2 p 2 i 4 ) } / (n 5 (p - i) 5 (p + i) 5 ) ; i * p 

2 2 ^ 

- 71 D - 3 . 

= ;i = p , 

■a a ’ r » 


M f P ° = Jcos(iTC^) sin(prt^)d^ = M 2 ®, 

0 

1 

M 3 p = J^cos (in^) sin (pJt^)d^ = M 2 ) . 
0 

M ip 2 = jVcos(iJt^)sin(pJt£)d£ = M 22 , 
0 


(B.153) 


(B.154) 


6) (B.155) 


(B. 156) 


(B.157) 


(B.158) 
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M” = jVcos (iit£) sin (pn£) d£ = M 23 , 
o 


M 34 = jVcos (iJi£) sin (pJt^)d^ = M 24 , 
o 
1 

M 4 ^ 1 = Jcos (i7i^) cos (pn^)d^ = 0; i *p 
o 

1 . 

= -;» = P. 


M 41 = f^cos (in!;) cos (pJt£) d^ = 

J Jt (p-i 


(p 2 + i 2 )((-l) p+ '-l) 

2 /_ . :\ 2 ’^ P 


(p-i) (p + i)" 


= -; 1=P 


M 42 = jVcos (in£) cos (pn£) d^ = 2 ^ P 

J n (p-i 


2 (p 2 + i 2 ) (-D P+i . : 

2 /_ . :,2 


; i ^ p 


(p-i) (p + i)' 
2n 2 p 2 + 3 . 

= — Vt ;i = p ■ 

12tc p 


M 43 = jVcos (iJt£) cos (pTt^)d^ 


3 [ft 2 (i 6 _ i 4 p 2 ~ i 2 P 4 + P 6 ) (-l) P + i -2(i 4 + 6iV + p 4 )((-l) p + , -l)] , ; ^ i 

n 4 (p-i) 4 (p + i) 4 


2 2,o 

Tf + ? ;l = p 
87tV 


(B.159) 


(B.160) 


(B.161) 


(B.162) 


(B.163) 


(B.164) 
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M“-JV 


cos (i7c£,) cos (p7i^) 


_ 4 [n 2 (i 6 - i 4 p 2 - i 2 p 4 + p 6 ) - 6 (f* + 6i i p i + p‘ t ) ] (-1) 
^ 4 (p-i) 4 (p + i) 4 

2*V(*V + 5) -15 . 

— ;i = p , 

20n 4 p 4 


: 4 . ,,i_i . _4. 


, p +1 


-;i 5>tp 


and finally, 


m 45 = J*^ 5 cos (inJ;) cos (pJt^) dJ; 
o 

= 5[n 4 (i 10 -3iV + 2i 6 p 4 + 2i 4 p 6 -3i 2 p 8 + p 10 ) (-l) P+i 

-12n 2 (i 8 + 4i 6 p 2 - 10i 4 p 4 + 4i 2 p 6 + p 8 ) (— l) p+l 

+ 24(i 6 +15i 4 p 2 +15i 2 p 4 + p 6 ) ((-l) p+i -l) ]/(Ji 6 (p-i) 6 (p + i) 6 );i^p 

JC 2 p 2 (2K 2 p 2+ 15) -45 . 

_ ’ 1 ~ P ’ 

247ip 

For clamped boundary conditions on the ends of the plate, and <J) p are given by 

= cos [ (i - 1) rc£] - cos [ (i + 1) 7t£] 

<t> P = COS [ (P- 1) ic£] “COS [ (P + l)Jl£] 


Making this substitution, the integrands appearing in eqs. (B.50), (B.85), and (B.124) are given by 


J 


= e (M 4 * 1 — M? 21 — M 431 + M 441 ) + (M 4 * 0 - M 420 — M 430 + M? 40 ) , 

3 T i T p ^ v ip ip ip ip 7 v ip ip ip ip 7 ’ 

b 


J 


^pd^-^Kp+l) (M 2 ; 2 -M 2 p 32 ) - (P -1) (M 322 - M 342 ) ] 
b 


-2ne[(p+l) (M 3 p ! - M 331 ) - (p — 1 ) (M 321 — M 341 ) ] 

-* I (P + 1) (M 3 p'° - M 330 ) - (p - 1) (M 320 - M'p 40 ) ] 

, '>^2 /x . 41! x .421 ..431 , X .44K , - , x .410 x .420 x .430 , x .440. 

+ 2e (M. -M. -M. +M )+2£(m -M. -M +M )- 

v ip Ip ip ip y v lp lp ip ip ' 


(B.165) 


(B.166) 


(B.167) 


(B.168) 


(B.169) 
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J l ^"d§ = -4*e 3 [(p+l) (M^ 2 -M” 2 ) - (p-1) (M 222 -M 2 p 42 )] 

0 

-Site 2 [ (p + 1) (Mf; - M 3 p 31 ) - (P - 1 ) (M 321 - M™) ] 

-4 it [ (p + 1) (M 2 ;°-M 2 p 30 ) - (P-1) (Mf^-Mfp 40 )] 

-*V [ (p + 1) : 2 (M 4 ; 3 - M 4 p 33 ) - (p - 1) 2 (M 4 ? - M“ 3 ) ] 

-3itV [ ( P + 1 ) 2 (M 4 p 12 - M 4 p 32 ) - (p - 1 ) 2 (M 4 p 22 - M£ 2 ) ] 

-e[3it 2 (p + l) 2 -2e 2 ] (M 4 p u -M 4 p 31 ) + e[3it 2 ( P - l) 2 -2e 2 ] (Mj'-M* 1 ) 
- [it 2 (p + 1 ) 2 — 2e 2 ] (M 4 ;°-M 4 p 30 ) + [it 2 (p- 1) 2 -2e 2 ] (M^-Mj 0 ) > 


f i^;<t, p d^=-ite 2 [(i + i) (m 2 ; 2 -m 2 p 22 ) - (i-D (m 2 p 32 -m 2 ; 2 )] 

o b 

-2ne[ (i+ 1) (M^’-M 2 ; 1 ) -(1-1) (M 2 p 3, -M 2 ;‘)] (B.171) 

-tc [ (i + 1) (M 2 p ‘ - - M 2 p 20 ) - (i - 1 ) (M™ - M* 0 ) ] 

+ 2e 2 ( M 4 p u - M 421 - M 4 p 3] + Mj* ) + 2e ( M 4 ” - M«° - M 4 p 30 + M™) - 


d£ = itV [ (i + 1) ( ( P + 1) Mjp 3 - (P - 1) Mj 123 ) - (i - 1) ( (P + 1) M'p 33 - (p - 1)M’ 43 ) ] 
o b2 

+ 3itV[(i+l) ((p+DMjp 2 - (p- 1) M- 22 ) - (i - 1) ((p+l)M* p 32 - (p- 1) mJ 42 ) ] 

+ 3ltV((i+l) ((p+ (p- 1)M‘ p 21 ) - (i-1) ( (p+ 1)M' 31 - (p- 1)M? 41 )] 

+ it 2 [(i + 1) ((p+l)M?p°- (p-1) m' 20 ) - (i-1) ((p+lJM'p 30 - (p-IJM.'p 40 )] 

-2 lie 3 [ (i + 1 ) (M 2 ; 2 - M 222 ) - (i - 1) (M 232 - M 2 ; 2 ) ] (B 

-4ite 2 [ (i + 1) (M 2 ; 1 -M 221 ) - (i - 1) (M 231 - M 241 ) ] 


-2ite [ (i + 1) (M 2 ;°-M 2 p 20 ) - (i-1) (M 230 -M 2 ^)] 

-2ite 3 [ (p + 1) (M 3 p 12 - M 332 ) - (p - 1) (M 3 p 22 - M 3 ; 2 ) ] 
-4 ice 2 [ (p + 1 ) (M 3 ; 1 - M 3 p 31 ) - (p - 1 ) (M 3 p 21 - M 341 ) ] 


-2ite [ (p + 1) (M 3 p ° - M 3 p 30 ) - ( p - 1 ) (M 3 p 20 - M*°) ] 

+ 4e 3 (M 4 p ' - M 421 - M 4 p 31 + ) + 4E 2 ( M 4 ; 0 - M 420 - M 4 * + M*°) • 
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f = 4 n 2 e 4 [ (i + 1 ) ((p+l)M^ 3 - (p-l)M‘ p 23 ) - (i-1) ((p + 1)M? 33 - (p - 1 ) m' 43 ) ] 

J o b 

+ 12 JtV [ (i + 1) ( (p + 1) M^ 12 - (p - 1) Mj 122 ) - (i - 1) ( (p + 1 ) M’ 32 - (p - 1) M, 142 ) ] 

+ 12JtV[(i+l) ( (p + 1) M ip U - (p - 1) Mj 121 ) - (i-1) ((p+DMj 1 - (p-DMj 1 )] 

+ 4n 2 e [ (i + 1) ( (p + 1) M” 0 - (p - 1 ) M™) - (i - 1 ) ( ( P + 1 ) Ml 30 - (p - 1) m" 0 ) ] 

+ rtV [ (p + 1) 2 ( (i + 1) M 2 p 4 - (i - 1 ) M 234 ) - (p - l) 2 ( (i + 1) M 224 - (i - 1) M 2 ? 4 ) ] 

+ 4jt¥ [ (p + 1) 2 ( (i + 1) M 213 - (i - 1) M 233 ) - (p - 1) 2 ( (i + 1 ) M 223 - (i - 1) M 243 ) ] 

+ 6ji 3 £ 2 [ (p + 1) 2 ( (i + 1) M 2 , 12 - (i - 1) M 232 ) - (p - 1) 2 ( (i + 1 ) M 222 - (i - 1) M 242 ) ] 

+ 47i 3 e [ (p + 1) 2 ( (i + l)M? p " - (i - 1) M 231 ) - (p - l) 2 ( (i + I) M 221 - (i - 1) M 241 ) ] 

+ 7T 3 [ ( p + 1 ) 2 ( ( i + 1 ) M 2 p ' 0 - ( i - 1 ) Mfp 30 ) - ( P - 1 ) 2 ( ( i + 1 ) M 2 p 20 - ( i - 1 ) M fp 40 ) ] 

— 2tc£ 4 [ (i + 1) (M 2 p 12 - M 222 ) - (i - 1) (M^ 2 - M 242 ) ] f®- 173 ) 

— 47t£ 3 [ (i + 1) (M 211 - M 221 ) - (i - 1) (M 231 - M 241 ) ] 

-2JI£ 2 [ (i + 1 ) (M™ - M 230 ) - (i - 1 ) (M 220 - Mjp 40 ) ] 

-8JtE 4 [ (p+ 1) (M 3 p 12 -M 332 ) - (p- 1) (M 322 - M 342 ) ] 

-167t£ 3 [ (p+ 1) (Mf“ -M 331 ) - (p- 1) (M 321 -M 341 ) ] 

-8tc£ 2 [ (p + 1) (M 3 p '° - M 330 ) - (p - 1) (M - Mj 0 ) ] 

-2ti 2 £ 4 [ (p + 1) 2 (M 4 ; 3 - M 433 ) - (p - 1) 2 (M 4 * - M^ 3 ) ] 

-67tV [ (p + 1) 2 (M 4 p 12 - M 432 ) - (p - 1) 2 (M 422 - M 4 ^ 2 ) ] 

-6* V [ ( p + 1 ) 2 (M 41 ' 1 - M 4 p 31 ) - (p - 1 ) : 2 (M 421 - M*' ) ] 

-2^ 2 £ [ (p + 1 ) 2 (M 4 p 10 - M 4 p 30 ) - (p - 1 ) : 2 (M 420 - M J°) ] 

J = — 4it£ 3 [ (i + 1) (M 2 ; 2 - M 222 ) - (i - 1 ) (M 232 - M 242 ) ] 

0 

-87t£ 2 [ (i+ 1) (M 2 p n - M 221 ) - (i - 1) (M 231 - M 241 ) ] 

— 47t£[ (i+ 1) (M 210 -M 220 ) - (i- 1) (M 230 -M 240 )] 

v * v ip ip ’ v } v ip ip ' J (B.174) 

-JtV [ ( i + 1 ) 2 (M 413 - M 423 ) - (i - 1 ) 2 (M 4 p 33 - M^ 3 ) ] 

"3*V ( ( i + 1 ) 2 (Mj 2 - M 422 ) - (i - 1 ) ' 2 (M 4 p 32 - M“ 2 ) ] 

- £ [3ti 2 (i + 1) 2 - 2 e 2 ] (M 4 ” - M 421 ) + e [37t 2 (i - 1 ) ' 2 - 2 e 2 ] (M 431 - m£‘) 

- [7t 2 (i+l) 2 -2£ 2 ] (M 4 p 10 -M 420 ) + [jc 2 (i- 1) 2 -2 e 2 ] 
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= 47 t 2 e 4 [ (p + 1) ( (i + 1) M jp 3 - (i - 1 ) M, 133 ) - (p - 1 ) ( (i + 1 ) M, 123 - (i - 1 ) M; p 43 ) ] 

b 

o 

+ 12jcV [ (p + 1 ) ( (i + 1) M” 2 - (i - 1) m' 32 ) - (p - 1) ( (i + 1 ) M’ p 22 - (i - 1 ) m' 42 ) ] 

+ 12rc 2 £ 3 [ (p + 1) ( (i + 1) M-p 1 - (i - 1) M? 31 ) - (p - 1) ( (i + 1)M- 21 - (• - 1) M* 41 ) ] 

+ 4ji 2 e [ (p+ 1) ( (i+ 1) M'p 10 - (i-l)M™) -(p-1) ((i + DM, 120 - ( i - 1 ) Mf) ] 

-Sue 4 [ (i + 1) (M 2 p 12 - Mjf) - (i - 1) (M 232 - M 2 p 42 ) ] 

-16K£ 3 [(i + l) (M 2 ;‘-M 2 p 21 ) - (i-1) (Mj'-Mj 1 )] 

-87CE 2 [ (i + 1 ) (M 2 ; 0 - MjJ°) - (i - 1 ) (M*° - M*°) ] 

+ Jt 3 £ 4 [( P+ l)((i + l) 2 M 3 i ; 4 -(i-l) 2 Mf p 34 )-(p-l) ( (i + 1)2 M 3 p 24 -(i-l) 2 Mf p 44 )] 

+ 4K 3 e 3 [(p + l)((i + l) 2 M 3 p ,3 -(i-l) 2 M 333 )-( P -l)((i + l) 2 M 3 p 23 -(i-l) 2 M 3 p 43 )] 

+ 67t 3 e 2 [ (p+ 1) ( (i+ l) 2 M 3 p 12 - (i - 1) 2 M 332 ) - (p-1) ((i+l) 2 M 322 - (i - 1) 2 M 342 ) ] (B ‘ 175) 
+ 4 tc 3 £ [ (p + 1 ) ( ( i + 1 ) 'M 3 ! 1 1 - ( i - 1 ) 2 M 331 ) - ( p - 1 ) ( ( i + 1 ) 2 M 32 ‘ - (i - 1 ) 2 M*‘ ) ] 

+ k 3 [ ( P + 1) ( (i + 1) : 2 m 3 ;° - a - 1) 2 m 3 p 30 ) - ( P - 1) ( (i + d : 2 m ™ - (i - ir 2 m 3 p 40 ) ] 

-2*8 4 ((p+l) (M 3 p ,2 -M 332 ) - (p-1) (M 3 p 22 -M 3 p 42 )] 

-47t£ 3 [ (p + 1) (M 3 ” - M 331 ) - (p - 1) (M 3 p 21 - M 341 ) ] 

2jc£ 2 [ (P+1) (M™ - M 330 ) - (p - 1) (M* 0 - Mfp 40 ) ] 

-271 2 e 4 [ (i + 1 ) 2 (M 4 ; 3 - M 4 f) - (i - 1) : 2 (M 4 p 33 - M 4 ^ 3 ) ] 

-6rcV [ (i + 1 ) 2 (M 4 ; 2 - M 422 ) - (i - 1 ) 2 (M 4 p 32 - M 4 ^ 2 ) ] 

-67t 2 £ 2 [ (i + 1) ' 2 (M 4 ; 1 - M 4 p 21 ) - (i - 1) 2 (M 4 p 31 - M* 1 ) ] 

-2n 2 £ [ (i + 1 ) 2 (m 4 ; 0 - m; 20 ) - a - 1 ) 2 (m 4 p 30 - mj°) j • 


The M ip are the integrals over x, which have similar form to Eqs. B.145 - B.166. 
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Appendix C 

Details of Assembling Stiffness Matrices 


The total and geometric stiffness matrices, [K T ] and [K G ] , are each composed of four smaller 
submatrices. Specifically, 


[K t ] 


, EE 1 , OE 
k T i k x 

w EOI w°° 

l\y | Iv'j' 


[k g ] = 


* EE , OE 
k G J_ k _G_ 

. EOt .OO 
K r. » K r, 


(C.i) 


Each element of a submatrix is the result of integrating a combination of derivatives of w, in the series 
approximation for the out-of-plane displacement of the plate, over the surface area of the plate. As such, 
each element is a function of the plate geometry, material properties, and the coefficients i, j, p, and q. The i 
and p indices are from the sets of functions of x, the coordinate along the length of the plate 
( < ^ >i I °° 12 > j 2 ) f° r w an( ^ 8w, respectively. The j and q indices are from the sets of functions of y, 

the coordinate across the width of the plate Vj (fj) |°° 1 2 > ^ (gj) |°°_ ^ 2 , 

V q (f q ) | °° { 2 »r q (g q ) |°° i ^ . Since <j>, \\f and T are trigonometric functions, the resulting integrations 

over the period of these functions yield different results that depend on the particular combination of (i,p) 
and (j,q). Special cases occur because of the orthogonality of the trigonometric functions, and these special 
cases are referred to as degenerate cases. As examples, consider the following cases for simply supported 
ends: 


and 


1 


cos (irc£) cos (pTt^)d^ 


-r. 


0 


i = p 
i*p 


(C.2) 
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b(4> 


J r|cos (2j - 


1) 


rcri 

2btfj. 


sin 


(2q-l) 


Ttri 

2b (£,) 


dr| 


-b ($) 


b 2 

(jt(2j-l)) 

~b 2 (2j — 1) (— l) j + q 
- 7C (q — j) (q+j- 1) 


j = q 

(C.3) 

j*q 


Thus, for simply supported cases there are four possibilities for combinations of i, j, p, and q: i = p and 
j s q,i = p and j ^ q , i & p and j = q , and finally, i ^ p and j ^ q . To demonstrate how this is handled in 

the computer code, the [K^ E ] {w E } multiplication is expanded. Remember that in each row of the matrix, 

p and q are fixed indices while i and j are free indices. The expression is 


[ K SJ = S S { A dpqVjq + B dP« 5 'P + C ‘ipq 5 iq + D .JPq } 

i=lj=l 

3 3 3 3 

= a„w' + X B ,„ W „ + S c . w *i, + E5X>. W » 

j = 1 i = 1 i = lj = 1 

i * o i * n i * n i 


(C.4) 


The Apq are coefficients resulting from integrations performed for i = p and j - q , the B p j q are coefficients 
resulting from integrations performed for i = p and j ^ q , the Cjpq are coefficients resulting from 
integrations performed for i * p and j = q , the D ijpq are coefficients resulting from integrations performed 
for i * p, and j * q . The submatrices are constructed row-by-row in the computer program with p and q 
being fixed values for each particular row. Four procedures are implemented to fill the columns of each row. 
The first process computes the single term Apq, which is always on the diagonal. The second process 
computes the terms B p j q . To do this, a loop over j, j = 1, 2, ..., N, is initiated (N being the number of terms 
in the series approximation for w). For each execution of this loop, as long as j * q, B pjq is computed and 
loaded into the appropriate column. This appropriate column is computed according to: 


column = row- q+ 1 first pass through loop 
column = column + 1 indexed for subsequent passes. 


The third process computes the terms C ipq . To do this, a loop over i, i = 1, 2, .... N, is initiated. For each 
execution of this loop, as long as i * p, is computed and loaded into the appropriate column. This 
appropriate column is computed according to: 
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column = q first pass through loop 

column = column + N indexed for subsequent passes. 


(C.6) 


The fourth process computes the terms D ip j q . To do this, two nested loops over i and j, i = 1, 2, ...» N and, 
j = 1 , 2, . . N are initiated. For each execution of these loops, as long as i * p and j * q , D ip j q is computed 
and loaded into the appropriate column. This appropriate column is computed according to 


column = 1 
column = column + 1 


first pass through loop 
indexed for subsequent passes. 


(C.7) 


For N=3, a completely filled submatrix having size N xN has the structure of equation (C.8). 







i = 1 

i = 1 

i = 1 

i = 2 

i = 2 

i = 2 

i = 3 

i = 3 

i = 3 






j-1 

j = 2 

j = 3 

3 = 1 

j = 2 

j = 3 

j-1 

j = 2 

J-3 

p = 

1 

q 

= 

l 

A n 

B 121 

B 131 

r 

^211 

^2121 

^2131 

^311 

D 3121 

^3131 

p = 

1 

q 

= 

2 

® 112 

A]2 

B 132 

^2112 

c 

^212 

^2131 

^3112 

C 3.2 

^3132 

p = 

1 

q 

= 

3 

B 113 

B 123 

Aj3 

^2113 

^2123 

C 213 

D 3113 

D 3.23 

c 

^313 

p - 

2 

q 

= 

1 

C 121 

^1221 

^1231 

A 2 i 

B 221 

B 231 

c 

v "321 

^3221 

^3231 

p = 

2 

q 

* 

2 

^1212 

c 

^122 

^1232 

B 212 

A 22 

B 232 

^3212 

c 

'"'322 

^3232 

p = 

2 

q 

= 

3 

^1213 

^1223 

c 

'“'123 

B 213 

B 223 

A 23 

D 3213 

D 3223 

^323 

p = 

3 

q 


1 

^131 

^1321 

^1331 

^231 

^2321 

D 2331 

A 31 

B 32, 

®33I 

p “ 

3 

q 

= 

2 

D 1312 

c 

^132 

^1332 

^2312 

c 

^232 

D 2 33 2 

B 312 

^32 

B 332 

p = 

3 

q 

= 

3 

Pl313 

^1323 

c 

'“'133 

D 2313 

^2323 

c 

^233 

B 313 

B 323 

A 33_ 


(C.8) 
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